24 #include "absl/container/flat_hash_map.h"
25 #include "absl/numeric/int128.h"
33 #include "ortools/glop/parameters.pb.h"
52 for (
const glop::ColIndex
col : non_zeros_) {
53 dense_vector_[
col] = IntegerValue(0);
55 dense_vector_.
resize(size, IntegerValue(0));
57 dense_vector_.
assign(size, IntegerValue(0));
59 for (
const glop::ColIndex
col : non_zeros_) {
60 is_zeros_[
col] =
true;
62 is_zeros_.
resize(size,
true);
70 dense_vector_[
col] = IntegerValue(add);
71 if (is_sparse_ && is_zeros_[
col]) {
72 is_zeros_[
col] =
false;
73 non_zeros_.push_back(
col);
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);
88 if (!
AddProductTo(multiplier, term.second, &dense_vector_[term.first])) {
92 if (
static_cast<double>(non_zeros_.size()) < threshold) {
97 for (
const std::pair<glop::ColIndex, IntegerValue> term : terms) {
98 if (!
AddProductTo(multiplier, term.second, &dense_vector_[term.first])) {
107 const std::vector<IntegerVariable>& integer_variables,
109 result->
vars.clear();
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);
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);
129 result->
ub = upper_bound;
132 std::vector<std::pair<glop::ColIndex, IntegerValue>>
134 std::vector<std::pair<glop::ColIndex, IntegerValue>> result;
136 std::sort(non_zeros_.begin(), non_zeros_.end());
137 for (
const glop::ColIndex
col : non_zeros_) {
138 const IntegerValue coeff = dense_vector_[
col];
142 const int size = dense_vector_.
size();
143 for (glop::ColIndex
col(0);
col < size; ++
col) {
144 const IntegerValue coeff = dense_vector_[
col];
154 : constraint_manager_(
model),
155 sat_parameters_(*(
model->GetOrCreate<SatParameters>())),
162 implied_bounds_processor_({}, integer_trail_,
165 expanded_lp_solution_(
171 if (sat_parameters_.use_branching_in_lp() ||
172 sat_parameters_.search_branching() == SatParameters::LP_SEARCH) {
173 compute_reduced_cost_averages_ =
true;
177 integer_trail_->RegisterReversibleClass(&rc_rev_int_repository_);
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;
186 << num_solves_by_status_[i];
192 DCHECK(!lp_constraint_is_registered_);
193 constraint_manager_.
Add(
ct);
200 for (
const IntegerVariable
var :
ct.vars) {
205 glop::ColIndex LinearProgrammingConstraint::GetOrCreateMirrorVariable(
206 IntegerVariable 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());
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;
220 if (
index >= expanded_lp_solution_.
size()) {
229 IntegerValue coeff) {
230 CHECK(!lp_constraint_is_registered_);
231 objective_is_defined_ =
true;
233 if (ivar != pos_var) coeff = -coeff;
236 const glop::ColIndex
col = GetOrCreateMirrorVariable(pos_var);
237 integer_objective_.push_back({
col, coeff});
238 objective_infinity_norm_ =
255 bool LinearProgrammingConstraint::CreateLpFromConstraintManager() {
258 infinity_norms_.
clear();
259 const auto& all_constraints = constraint_manager_.
AllConstraints();
263 integer_lp_.
push_back(LinearConstraintInternal());
264 LinearConstraintInternal& new_ct = integer_lp_.
back();
267 const int size =
ct.vars.size();
268 IntegerValue infinity_norm(0);
270 VLOG(1) <<
"Trivial infeasible bound in an LP constraint";
279 for (
int i = 0; i < size; ++i) {
281 IntegerVariable
var =
ct.vars[i];
282 IntegerValue coeff =
ct.coeffs[i];
288 new_ct.terms.push_back({GetOrCreateMirrorVariable(
var), coeff});
290 infinity_norms_.
push_back(infinity_norm);
293 std::sort(new_ct.terms.begin(), new_ct.terms.end());
298 for (
int i = 0; i < integer_variables_.size(); ++i) {
305 objective_infinity_norm_ = 0;
306 for (
const auto entry : integer_objective_) {
307 const IntegerVariable
var = integer_variables_[entry.first.value()];
309 integer_objective_offset_ +=
313 objective_infinity_norm_ =
315 integer_objective_[new_size++] = entry;
318 objective_infinity_norm_ =
320 integer_objective_.resize(new_size);
323 for (
const LinearConstraintInternal&
ct : integer_lp_) {
326 for (
const auto& term :
ct.terms) {
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];
348 glop::GlopParameters params;
349 params.set_cost_scaling(glop::GlopParameters::MEAN_COST_SCALING);
350 scaler_.
Scale(params, &lp_data_);
351 UpdateBoundsOfLpVariables();
356 if (sat_parameters_.polish_lp_solution()) {
358 for (
int i = 0; i < num_vars; ++i) {
359 const IntegerVariable cp_var = integer_variables_[i];
362 if (lb != 0 || ub != 1)
continue;
373 <<
" Managed constraints.";
377 LPSolveInfo LinearProgrammingConstraint::SolveLpForBranching() {
379 glop::BasisState basis_state = simplex_.
GetState();
381 const glop::Status status = simplex_.
Solve(lp_data_, time_limit_);
385 VLOG(1) <<
"The LP solver encountered an error: " << status.error_message();
394 info.new_obj_bound = IntegerValue(
395 static_cast<int64>(std::ceil(info.lp_objective - kCpEpsilon)));
400 void LinearProgrammingConstraint::FillReducedCostReasonIn(
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(
410 }
else if (rc < -kLpEpsilon) {
411 integer_reason->push_back(
419 bool LinearProgrammingConstraint::BranchOnVar(IntegerVariable positive_var) {
421 DCHECK(lp_solution_is_set_);
423 DCHECK_GT(std::abs(current_value - std::round(current_value)), kCpEpsilon);
426 integer_reason_.clear();
428 bool deductions_were_made =
false;
430 UpdateBoundsOfLpVariables();
432 const IntegerValue current_obj_lb = integer_trail_->
LowerBound(objective_cp_);
436 const glop::ColIndex lp_var = GetOrCreateMirrorVariable(positive_var);
440 if (current_value < current_lb || current_value > current_ub) {
445 const double new_ub = std::floor(current_value);
448 LPSolveInfo lower_branch_info = SolveLpForBranching();
458 positive_var, IntegerValue(std::ceil(current_value)));
459 if (!integer_trail_->
Enqueue(deduction, {}, integer_reason_)) {
462 deductions_were_made =
true;
463 }
else if (lower_branch_info.new_obj_bound <= current_obj_lb) {
468 const double new_lb = std::ceil(current_value);
471 LPSolveInfo upper_branch_info = SolveLpForBranching();
475 return deductions_were_made;
482 positive_var, IntegerValue(std::floor(current_value)));
483 if (!integer_trail_->
Enqueue(deduction, {}, integer_reason_)) {
484 return deductions_were_made;
486 deductions_were_made =
true;
488 }
else if (upper_branch_info.new_obj_bound <= current_obj_lb) {
489 return deductions_were_made;
498 approximate_obj_lb = upper_branch_info.new_obj_bound;
500 approximate_obj_lb = lower_branch_info.new_obj_bound;
502 approximate_obj_lb =
std::min(lower_branch_info.new_obj_bound,
503 upper_branch_info.new_obj_bound);
508 if (approximate_obj_lb <= current_obj_lb)
return deductions_were_made;
511 const IntegerLiteral deduction =
513 if (!integer_trail_->
Enqueue(deduction, {}, integer_reason_)) {
514 return deductions_were_made;
521 DCHECK(!lp_constraint_is_registered_);
522 lp_constraint_is_registered_ =
true;
527 std::sort(integer_objective_.begin(), integer_objective_.end());
530 if (!sat_parameters_.add_lp_constraints_lazily()) {
533 if (!CreateLpFromConstraintManager()) {
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++) {
544 if (objective_is_defined_) {
557 optimal_constraints_.resize(rev_optimal_constraints_size_);
558 if (lp_solution_is_set_ && level < lp_solution_level_) {
559 lp_solution_is_set_ =
false;
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])] =
580 for (
const IntegerVariable
var : generator.
vars) {
583 cut_generators_.push_back(std::move(generator));
587 const std::vector<int>& watch_indices) {
588 if (!lp_solution_is_set_)
return Propagate();
598 for (
const int index : watch_indices) {
604 if (value < lb - kCpEpsilon || value > ub + kCpEpsilon)
return Propagate();
619 glop::ColIndex
var) {
624 IntegerVariable variable)
const {
625 return lp_solution_[
gtl::FindOrDie(mirror_lp_variable_, variable).value()];
629 IntegerVariable variable)
const {
630 return lp_reduced_cost_[
gtl::FindOrDie(mirror_lp_variable_, variable)
634 void LinearProgrammingConstraint::UpdateBoundsOfLpVariables() {
635 const int num_vars = integer_variables_.size();
636 for (
int i = 0; i < num_vars; i++) {
637 const IntegerVariable cp_var = integer_variables_[i];
645 bool LinearProgrammingConstraint::SolveLp() {
647 lp_at_level_zero_is_final_ =
false;
650 const auto status = simplex_.
Solve(lp_data_, time_limit_);
653 VLOG(1) <<
"The LP solver encountered an error: " << status.error_message();
657 average_degeneracy_.
AddData(CalculateDegeneracy());
659 VLOG(2) <<
"High average degeneracy: "
664 if (status_as_int >= num_solves_by_status_.size()) {
665 num_solves_by_status_.resize(status_as_int + 1);
667 num_solves_by_status_[status_as_int]++;
674 lp_solution_is_set_ =
true;
676 const int num_vars = integer_variables_.size();
677 for (
int i = 0; i < num_vars; i++) {
679 GetVariableValueAtCpScale(glop::ColIndex(i));
680 lp_solution_[i] =
value;
681 expanded_lp_solution_[integer_variables_[i]] =
value;
685 if (lp_solution_level_ == 0) {
686 level_zero_lp_solution_ = lp_solution_;
692 bool LinearProgrammingConstraint::AddCutFromConstraints(
693 const std::string&
name,
694 const std::vector<std::pair<RowIndex, IntegerValue>>& integer_multipliers) {
705 if (!ComputeNewLinearConstraint(integer_multipliers, &tmp_scattered_vector_,
707 VLOG(1) <<
"Issue, overflow!";
726 if (std::abs(activity -
ToDouble(cut_.
ub)) / norm > 1e-4) {
727 VLOG(1) <<
"Cut not tight " << activity <<
" <= " <<
ToDouble(cut_.
ub);
737 const IntegerVariable first_new_var(expanded_lp_solution_.
size());
738 CHECK_EQ(first_new_var.value() % 2, 0);
740 LinearConstraint copy_in_debug;
742 copy_in_debug = cut_;
752 std::vector<ImpliedBoundsProcessor::SlackInfo> ib_slack_infos;
754 false, first_new_var,
755 expanded_lp_solution_, &cut_, &ib_slack_infos);
757 cut_, ib_slack_infos));
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) {
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);
777 tmp_lp_values_.push_back(expanded_lp_solution_[
var]);
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;
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);
801 const IntegerValue diff(
802 CapSub(integer_lp_[
row].ub.value(), integer_lp_[
row].lb.value()));
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);
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));
814 bool at_least_one_added =
false;
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());
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(
835 absl::StrCat(
"num_lifted_booleans=",
837 first_new_var, first_slack, ib_slack_infos, &cut_);
839 return at_least_one_added;
842 bool LinearProgrammingConstraint::PostprocessAndAddCut(
843 const std::string&
name,
const std::string& info,
844 IntegerVariable first_new_var, IntegerVariable first_slack,
845 const std::vector<ImpliedBoundsProcessor::SlackInfo>& ib_slack_infos,
846 LinearConstraint* cut) {
850 double activity = 0.0;
851 for (
int i = 0; i < cut->vars.size(); ++i) {
852 if (cut->vars[i] < first_new_var) {
854 ToDouble(cut->coeffs[i]) * expanded_lp_solution_[cut->vars[i]];
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);
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];
874 if (
var < first_new_var) {
875 const glop::ColIndex
col =
878 tmp_scattered_vector_.
Add(
col, cut->coeffs[i]);
880 tmp_scattered_vector_.
Add(
col, -cut->coeffs[i]);
886 if (
var < first_slack) {
887 const IntegerValue multiplier = cut->coeffs[i];
888 const int index = (
var.value() - first_new_var.value()) / 2;
891 std::vector<std::pair<ColIndex, IntegerValue>> terms;
892 for (
const std::pair<IntegerVariable, IntegerValue>& term :
893 ib_slack_infos[
index].terms) {
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];
917 multiplier, integer_lp_[
row].terms)) {
923 if (!
AddProductTo(multiplier, tmp_slack_bounds_[slack_index], &cut_ub)) {
930 VLOG(1) <<
"Overflow in slack removal.";
934 VLOG(3) <<
" num_slack: " << num_slack;
940 const std::string extra_info =
941 absl::StrCat(info,
" num_ib_substitutions=", ib_slack_infos.size());
943 const double new_violation =
945 if (std::abs(violation - new_violation) >= 1e-4) {
946 VLOG(1) <<
"Violation discrepancy after slack removal. "
947 <<
" before = " << violation <<
" after = " << new_violation;
951 return constraint_manager_.
AddCut(*cut,
name, expanded_lp_solution_,
959 void LinearProgrammingConstraint::AddCGCuts() {
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) {
965 const Fractional lp_value = GetVariableValueAtCpScale(basis_col);
973 if (std::abs(lp_value - std::round(lp_value)) < 0.01)
continue;
977 if (basis_col >= integer_variables_.size())
continue;
982 double magnitude = 0.0;
983 lp_multipliers.clear();
985 if (lambda.non_zeros.empty()) {
986 for (RowIndex
row(0);
row < num_rows; ++
row) {
988 if (std::abs(
value) < kZeroTolerance)
continue;
994 VLOG(1) <<
"BASIC row not expected! " <<
value;
999 lp_multipliers.push_back({
row,
value});
1002 for (
const ColIndex
col : lambda.non_zeros) {
1004 const double value = lambda.values[
col];
1005 if (std::abs(
value) < kZeroTolerance)
continue;
1009 VLOG(1) <<
"BASIC row not expected! " <<
value;
1014 lp_multipliers.push_back({
row,
value});
1017 if (lp_multipliers.empty())
continue;
1020 for (
int i = 0; i < 2; ++i) {
1026 for (std::pair<RowIndex, double>& p : lp_multipliers) {
1027 p.second = -p.second;
1033 integer_multipliers =
1034 ScaleLpMultiplier(
false,
1035 lp_multipliers, &scaling, 52);
1036 AddCutFromConstraints(
"CG", integer_multipliers);
1044 void RandomPick(
const std::vector<RowIndex>&
a,
const std::vector<RowIndex>&
b,
1045 ModelRandomGenerator* random,
1046 std::vector<std::pair<RowIndex, RowIndex>>* output) {
1047 if (
a.empty() ||
b.empty())
return;
1048 for (
const RowIndex
row :
a) {
1049 const RowIndex other =
b[absl::Uniform<int>(*random, 0,
b.size())];
1051 output->push_back({
row, other});
1056 template <
class ListOfTerms>
1057 IntegerValue GetCoeff(ColIndex
col,
const ListOfTerms& terms) {
1058 for (
const auto& term : terms) {
1059 if (term.first ==
col)
return term.second;
1061 return IntegerValue(0);
1066 void LinearProgrammingConstraint::AddMirCuts() {
1082 integer_variables_.size(), IntegerValue(0));
1083 SparseBitset<ColIndex> non_zeros_(ColIndex(integer_variables_.size()));
1088 std::vector<std::pair<RowIndex, IntegerValue>> base_rows;
1090 for (RowIndex
row(0);
row < num_rows; ++
row) {
1097 base_rows.push_back({
row, IntegerValue(1)});
1101 base_rows.push_back({
row, IntegerValue(-1)});
1122 std::vector<double> weights;
1124 std::vector<std::pair<RowIndex, IntegerValue>> integer_multipliers;
1125 for (
const std::pair<RowIndex, IntegerValue>& entry : base_rows) {
1135 integer_multipliers = {entry};
1136 if (AddCutFromConstraints(
"MIR_1", integer_multipliers)) {
1141 for (
const ColIndex
col : non_zeros_.PositionsSetAtLeastOnce()) {
1142 dense_cut[
col] = IntegerValue(0);
1144 non_zeros_.SparseClearAll();
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;
1156 used_rows.
assign(num_rows.value(),
false);
1157 used_rows[entry.first] =
true;
1162 const int kMaxAggregation = 5;
1163 for (
int i = 0; i < kMaxAggregation; ++i) {
1166 IntegerValue max_magnitude(0);
1168 std::vector<ColIndex> col_candidates;
1169 for (
const ColIndex
col : non_zeros_.PositionsSetAtLeastOnce()) {
1170 if (dense_cut[
col] == 0)
continue;
1173 const int col_degree =
1175 if (col_degree <= 1)
continue;
1180 const IntegerVariable
var = integer_variables_[
col.value()];
1181 const double lp_value = expanded_lp_solution_[
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);
1190 if (col_candidates.empty())
break;
1192 const ColIndex var_to_eliminate =
1193 col_candidates[std::discrete_distribution<>(weights.begin(),
1194 weights.end())(*random_)];
1197 std::vector<RowIndex> possible_rows;
1200 const RowIndex
row = entry.row();
1208 if (used_rows[
row])
continue;
1209 used_rows[
row] =
true;
1217 bool add_row =
false;
1220 if (entry.coefficient() > 0.0) {
1221 if (dense_cut[var_to_eliminate] < 0) add_row =
true;
1223 if (dense_cut[var_to_eliminate] > 0) add_row =
true;
1228 if (entry.coefficient() > 0.0) {
1229 if (dense_cut[var_to_eliminate] > 0) add_row =
true;
1231 if (dense_cut[var_to_eliminate] < 0) add_row =
true;
1236 weights.push_back(row_weights[
row]);
1239 if (possible_rows.empty())
break;
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);
1248 IntegerValue mult1 = -to_combine_coeff;
1249 IntegerValue mult2 = dense_cut[var_to_eliminate];
1256 const IntegerValue gcd = IntegerValue(
1266 for (std::pair<RowIndex, IntegerValue>& entry : integer_multipliers) {
1267 max_magnitude =
std::max(max_magnitude, entry.second);
1269 if (
CapAdd(
CapProd(max_magnitude.value(), std::abs(mult1.value())),
1271 std::abs(mult2.value()))) ==
kint64max) {
1275 for (std::pair<RowIndex, IntegerValue>& entry : integer_multipliers) {
1276 entry.second *= mult1;
1278 integer_multipliers.push_back({row_to_combine, mult2});
1281 if (AddCutFromConstraints(absl::StrCat(
"MIR_", i + 2),
1282 integer_multipliers)) {
1288 if (i + 1 == kMaxAggregation)
break;
1290 for (ColIndex
col : non_zeros_.PositionsSetAtLeastOnce()) {
1291 dense_cut[
col] *= mult1;
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;
1304 void LinearProgrammingConstraint::AddZeroHalfCuts() {
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]);
1319 for (glop::RowIndex
row(0);
row < integer_lp_.size(); ++
row) {
1327 row, integer_lp_[
row].terms, integer_lp_[
row].lb, integer_lp_[
row].ub);
1329 for (
const std::vector<std::pair<RowIndex, IntegerValue>>& multipliers :
1337 AddCutFromConstraints(
"ZERO_HALF", multipliers);
1341 void LinearProgrammingConstraint::UpdateSimplexIterationLimit(
1342 const int64 min_iter,
const int64 max_iter) {
1343 if (sat_parameters_.linearization_level() < 2)
return;
1344 const int64 num_degenerate_columns = CalculateDegeneracy();
1346 if (num_cols <= 0) {
1350 const int64 decrease_factor = (10 * num_degenerate_columns) / num_cols;
1355 if (is_degenerate_) {
1358 next_simplex_iter_ *= 2;
1361 if (is_degenerate_) {
1362 next_simplex_iter_ /=
std::max(
int64{1}, 2 * decrease_factor);
1366 next_simplex_iter_ = num_cols / 40;
1369 next_simplex_iter_ =
1374 UpdateBoundsOfLpVariables();
1379 if ( (
false) && objective_is_defined_) {
1388 static_cast<double>(integer_trail_->
UpperBound(objective_cp_).value() +
1389 100.0 * kCpEpsilon));
1398 parameters.set_max_number_of_iterations(2000);
1400 parameters.set_max_number_of_iterations(next_simplex_iter_);
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);
1410 if (!SolveLp())
return true;
1413 const int max_cuts_rounds =
1415 ? sat_parameters_.max_cut_rounds_at_level_zero()
1419 cuts_round < max_cuts_rounds) {
1423 if (!integer_lp_.empty()) {
1426 expanded_lp_solution_);
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();
1439 if (!cut_generators_.empty() &&
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_);
1448 expanded_lp_solution_, &constraint_manager_);
1452 if (constraint_manager_.
ChangeLp(expanded_lp_solution_, &state)) {
1454 if (!CreateLpFromConstraintManager()) {
1458 if (!SolveLp())
return true;
1460 VLOG(1) <<
"Relaxation improvement " << old_obj <<
" -> "
1467 lp_at_level_zero_is_final_ =
true;
1475 if (sat_parameters_.use_exact_lp_reason()) {
1476 if (!FillExactDualRayReason())
return true;
1485 UpdateSimplexIterationLimit(10, 1000);
1488 if (objective_is_defined_ &&
1495 const IntegerValue approximate_new_lb(
1496 static_cast<int64>(std::ceil(relaxed_optimal_objective - kCpEpsilon)));
1500 if (sat_parameters_.use_exact_lp_reason()) {
1501 if (!ExactLpReasonning())
return false;
1504 const IntegerValue propagated_lb =
1506 if (approximate_new_lb > propagated_lb) {
1507 VLOG(2) <<
"LP objective [ " <<
ToDouble(propagated_lb) <<
", "
1509 <<
" ] approx_lb += "
1510 <<
ToDouble(approximate_new_lb - propagated_lb) <<
" gap: "
1511 << integer_trail_->
UpperBound(objective_cp_) - propagated_lb;
1514 FillReducedCostReasonIn(simplex_.
GetReducedCosts(), &integer_reason_);
1515 const double objective_cp_ub =
1517 ReducedCostStrengtheningDeductions(objective_cp_ub -
1518 relaxed_optimal_objective);
1519 if (!deductions_.empty()) {
1520 deductions_reason_ = integer_reason_;
1521 deductions_reason_.push_back(
1526 if (approximate_new_lb > integer_trail_->
LowerBound(objective_cp_)) {
1529 if (!integer_trail_->
Enqueue(deduction, {}, integer_reason_)) {
1535 if (!deductions_.empty()) {
1536 const int trail_index_with_same_reason = integer_trail_->
Index();
1538 if (!integer_trail_->
Enqueue(deduction, {}, deductions_reason_,
1539 trail_index_with_same_reason)) {
1549 CHECK(lp_solution_is_set_);
1552 lp_solution_is_integer_ =
true;
1553 const int num_vars = integer_variables_.size();
1554 for (
int i = 0; i < num_vars; i++) {
1557 if (std::abs(lp_solution_[i] - std::round(lp_solution_[i])) >
1559 lp_solution_is_integer_ =
false;
1563 if (compute_reduced_cost_averages_) {
1564 UpdateAverageReducedCosts();
1568 if (sat_parameters_.use_branching_in_lp() && objective_is_defined_ &&
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_) {
1578 count_since_last_branching_ = 0;
1579 bool branching_successful =
false;
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];
1591 if (std::abs(current_value - std::round(current_value)) <= kCpEpsilon) {
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);
1611 branching_vars.insert(iterator, branching_var);
1612 if (branching_vars.size() > max_num_branches) {
1613 branching_vars.resize(max_num_branches);
1617 for (
const std::pair<double, IntegerVariable>& branching_var :
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;
1628 if (!branching_successful) {
1629 branching_frequency_ *= 2;
1638 IntegerValue LinearProgrammingConstraint::GetImpliedLowerBound(
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];
1653 bool LinearProgrammingConstraint::PossibleOverflow(
1654 const LinearConstraint& constraint) {
1655 IntegerValue lower_bound(0);
1656 const int size = constraint.vars.size();
1657 for (
int i = 0; i < size; ++i) {
1658 const IntegerVariable
var = constraint.vars[i];
1659 const IntegerValue coeff = constraint.coeffs[i];
1661 const IntegerValue
bound = coeff > 0
1668 const int64 slack =
CapAdd(lower_bound.value(), -constraint.ub.value());
1677 absl::int128 FloorRatio128(absl::int128 x, IntegerValue positive_div) {
1678 absl::int128 div128(positive_div.value());
1679 absl::int128 result = x / div128;
1680 if (result * div128 > x)
return result - 1;
1686 void LinearProgrammingConstraint::PreventOverflow(LinearConstraint* constraint,
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 =
1698 const double prod2 =
1703 const double max_value =
std::max({sum_max, -sum_min, sum_max - sum_min});
1705 const IntegerValue divisor(std::ceil(std::ldexp(max_value, -max_pow)));
1706 if (divisor <= 1)
return;
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);
1727 const absl::int128 remainder =
1728 absl::int128(old_coeff.value()) -
1729 absl::int128(new_coeff.value()) * absl::int128(divisor.value());
1735 if (new_coeff == 0)
continue;
1736 constraint->vars[new_size] = constraint->vars[i];
1737 constraint->coeffs[new_size] = new_coeff;
1740 constraint->vars.resize(new_size);
1741 constraint->coeffs.resize(new_size);
1743 constraint->ub = IntegerValue(
static_cast<int64>(
1744 FloorRatio128(absl::int128(constraint->ub.value()) - adjust, divisor)));
1749 void LinearProgrammingConstraint::SetImpliedLowerBoundReason(
1750 const LinearConstraint& terms, IntegerValue slack) {
1751 integer_reason_.clear();
1752 std::vector<IntegerValue> magnitudes;
1753 const int size = terms.vars.size();
1754 for (
int i = 0; i < size; ++i) {
1755 const IntegerVariable
var = terms.vars[i];
1756 const IntegerValue coeff = terms.coeffs[i];
1759 magnitudes.push_back(coeff);
1762 magnitudes.push_back(-coeff);
1773 std::vector<std::pair<RowIndex, IntegerValue>>
1774 LinearProgrammingConstraint::ScaleLpMultiplier(
1775 bool take_objective_into_account,
1776 const std::vector<std::pair<RowIndex, double>>& lp_multipliers,
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;
1786 if (std::abs(lp_multi) < kZeroTolerance)
continue;
1802 tmp_cp_multipliers_.push_back({
row, cp_multi});
1803 max_sum +=
ToDouble(infinity_norms_[
row]) * std::abs(cp_multi);
1808 if (take_objective_into_account) {
1809 max_sum +=
ToDouble(objective_infinity_norm_);
1813 std::vector<std::pair<RowIndex, IntegerValue>> integer_multipliers;
1814 if (max_sum == 0.0) {
1816 return integer_multipliers;
1821 const double threshold = std::ldexp(1, max_pow) / max_sum;
1822 if (threshold < 1.0) {
1825 return integer_multipliers;
1827 while (2 * *scaling <= threshold) *scaling *= 2;
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});
1836 return integer_multipliers;
1839 bool LinearProgrammingConstraint::ComputeNewLinearConstraint(
1840 const std::vector<std::pair<RowIndex, IntegerValue>>& integer_multipliers,
1841 ScatteredIntegerVector* scattered_vector, IntegerValue* upper_bound)
const {
1844 scattered_vector->ClearAndResize(integer_variables_.size());
1848 for (
const std::pair<RowIndex, IntegerValue> term : integer_multipliers) {
1849 const RowIndex
row = term.first;
1850 const IntegerValue multiplier = term.second;
1854 if (!scattered_vector->AddLinearExpressionMultiple(
1855 multiplier, integer_lp_[
row].terms)) {
1860 const IntegerValue
bound =
1861 multiplier > 0 ? integer_lp_[
row].ub : integer_lp_[
row].lb;
1869 void LinearProgrammingConstraint::AdjustNewLinearConstraint(
1870 std::vector<std::pair<glop::RowIndex, IntegerValue>>* integer_multipliers,
1871 ScatteredIntegerVector* scattered_vector, IntegerValue* upper_bound)
const {
1872 const IntegerValue kMaxWantedCoeff(1e18);
1873 for (std::pair<RowIndex, IntegerValue>& term : *integer_multipliers) {
1874 const RowIndex
row = term.first;
1875 const IntegerValue multiplier = term.second;
1876 if (multiplier == 0)
continue;
1880 IntegerValue negative_limit = kMaxWantedCoeff;
1881 IntegerValue positive_limit = kMaxWantedCoeff;
1885 if (integer_lp_[
row].ub != integer_lp_[
row].lb) {
1886 if (multiplier > 0) {
1887 negative_limit =
std::min(negative_limit, multiplier);
1889 positive_limit =
std::min(positive_limit, -multiplier);
1894 const IntegerValue row_bound =
1895 multiplier > 0 ? integer_lp_[
row].ub : integer_lp_[
row].lb;
1896 if (row_bound != 0) {
1900 const IntegerValue limit2 =
1902 if ((*upper_bound > 0) == (row_bound > 0)) {
1903 positive_limit =
std::min(positive_limit, limit1);
1904 negative_limit =
std::min(negative_limit, limit2);
1906 negative_limit =
std::min(negative_limit, limit1);
1907 positive_limit =
std::min(positive_limit, limit2);
1919 double positive_diff =
ToDouble(row_bound);
1920 double negative_diff =
ToDouble(row_bound);
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);
1931 const IntegerVariable
var = integer_variables_[
col.value()];
1938 const IntegerValue current = (*scattered_vector)[
col];
1940 const IntegerValue overflow_limit(
1942 positive_limit =
std::min(positive_limit, overflow_limit);
1943 negative_limit =
std::min(negative_limit, overflow_limit);
1960 const IntegerValue current_magnitude =
IntTypeAbs(current);
1961 const IntegerValue other_direction_limit =
FloorRatio(
1963 ? kMaxWantedCoeff +
std::min(current_magnitude,
1965 : current_magnitude,
1967 const IntegerValue same_direction_limit(
FloorRatio(
1968 std::max(IntegerValue(0), kMaxWantedCoeff - current_magnitude),
1970 if ((current > 0) == (coeff > 0)) {
1971 negative_limit =
std::min(negative_limit, other_direction_limit);
1972 positive_limit =
std::min(positive_limit, same_direction_limit);
1974 negative_limit =
std::min(negative_limit, same_direction_limit);
1975 positive_limit =
std::min(positive_limit, other_direction_limit);
1979 const IntegerValue implied = current > 0 ? lb : ub;
1990 IntegerValue to_add(0);
1991 if (positive_diff <= -1.0 && positive_limit > 0) {
1992 to_add = positive_limit;
1994 if (negative_diff >= 1.0 && negative_limit > 0) {
1997 std::abs(
ToDouble(negative_limit) * negative_diff) >
1998 std::abs(
ToDouble(positive_limit) * positive_diff)) {
1999 to_add = -negative_limit;
2003 term.second += to_add;
2004 *upper_bound += to_add * row_bound;
2008 CHECK(scattered_vector->AddLinearExpressionMultiple(
2009 to_add, integer_lp_[
row].terms));
2028 bool LinearProgrammingConstraint::ExactLpReasonning() {
2030 integer_reason_.clear();
2031 deductions_.clear();
2032 deductions_reason_.clear();
2038 std::vector<std::pair<RowIndex, double>> lp_multipliers;
2039 for (RowIndex
row(0);
row < num_rows; ++
row) {
2041 if (std::abs(
value) < kZeroTolerance)
continue;
2042 lp_multipliers.push_back({
row,
value});
2046 std::vector<std::pair<RowIndex, IntegerValue>> integer_multipliers =
2047 ScaleLpMultiplier(
true, lp_multipliers,
2051 if (!ComputeNewLinearConstraint(integer_multipliers, &tmp_scattered_vector_,
2053 VLOG(1) <<
"Issue while computing the exact LP reason. Aborting.";
2059 const IntegerValue obj_scale(std::round(scaling));
2060 if (obj_scale == 0) {
2061 VLOG(1) <<
"Overflow during exact LP reasoning. scaling=" << scaling;
2065 integer_objective_));
2067 AdjustNewLinearConstraint(&integer_multipliers, &tmp_scattered_vector_,
2072 LinearConstraint new_constraint;
2075 new_constraint.vars.push_back(objective_cp_);
2076 new_constraint.coeffs.push_back(-obj_scale);
2078 PreventOverflow(&new_constraint);
2079 DCHECK(!PossibleOverflow(new_constraint));
2082 IntegerSumLE* cp_constraint =
2083 new IntegerSumLE({}, new_constraint.vars, new_constraint.coeffs,
2084 new_constraint.ub, model_);
2088 optimal_constraints_.clear();
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();
2096 bool LinearProgrammingConstraint::FillExactDualRayReason() {
2099 std::vector<std::pair<RowIndex, double>> lp_multipliers;
2100 for (RowIndex
row(0);
row < ray.size(); ++
row) {
2102 if (std::abs(
value) < kZeroTolerance)
continue;
2103 lp_multipliers.push_back({
row,
value});
2105 std::vector<std::pair<RowIndex, IntegerValue>> integer_multipliers =
2106 ScaleLpMultiplier(
false, lp_multipliers,
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.";
2116 AdjustNewLinearConstraint(&integer_multipliers, &tmp_scattered_vector_,
2117 &new_constraint_ub);
2119 LinearConstraint new_constraint;
2121 integer_variables_, new_constraint_ub, &new_constraint);
2123 PreventOverflow(&new_constraint);
2124 DCHECK(!PossibleOverflow(new_constraint));
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;
2134 const IntegerValue slack = (implied_lb - new_constraint.ub) - 1;
2135 SetImpliedLowerBoundReason(new_constraint, slack);
2139 int64 LinearProgrammingConstraint::CalculateDegeneracy() {
2141 int num_non_basic_with_zero_rc = 0;
2142 for (glop::ColIndex i(0); i < num_vars; ++i) {
2144 if (rc != 0.0)
continue;
2148 num_non_basic_with_zero_rc++;
2151 is_degenerate_ = num_non_basic_with_zero_rc >= 0.3 * num_cols;
2152 return num_non_basic_with_zero_rc;
2155 void LinearProgrammingConstraint::ReducedCostStrengtheningDeductions(
2156 double cp_objective_delta) {
2157 deductions_.clear();
2162 const double lp_objective_delta =
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);
2171 if (rc == 0.0)
continue;
2172 const double lp_other_bound =
value + lp_objective_delta / rc;
2173 const double cp_other_bound =
2176 if (rc > kLpEpsilon) {
2178 const double new_ub = std::floor(cp_other_bound + kCpEpsilon);
2183 const IntegerValue new_ub_int(
static_cast<IntegerValue
>(new_ub));
2186 }
else if (rc < -kLpEpsilon) {
2188 const double new_lb = std::ceil(cp_other_bound - kCpEpsilon);
2190 const IntegerValue new_lb_int(
static_cast<IntegerValue
>(new_lb));
2191 deductions_.push_back(
2205 void AddOutgoingCut(
2206 int num_nodes,
int subset_size,
const std::vector<bool>& in_subset,
2207 const std::vector<int>& tails,
const std::vector<int>& heads,
2208 const std::vector<Literal>& literals,
2209 const std::vector<double>& literal_lp_values,
int64 rhs_lower_bound,
2211 LinearConstraintManager* manager, Model*
model) {
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;
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;
2241 if (num_optional_nodes_in + num_optional_nodes_out > 0) {
2243 rhs_lower_bound = 1;
2246 LinearConstraintBuilder outgoing(
model, IntegerValue(rhs_lower_bound),
2248 double sum_outgoing = 0.0;
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)));
2259 if (num_optional_nodes_in + num_optional_nodes_out > 0) {
2261 if (num_optional_nodes_in == subset_size &&
2262 (optional_loop_in == -1 ||
2263 literal_lp_values[optional_loop_in] > 1.0 - 1e-6)) {
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)) {
2273 if (num_optional_nodes_in == subset_size) {
2275 outgoing.AddLiteralTerm(literals[optional_loop_in], IntegerValue(1)));
2276 sum_outgoing += literal_lp_values[optional_loop_in];
2280 if (num_optional_nodes_out == num_nodes - subset_size) {
2281 CHECK(outgoing.AddLiteralTerm(literals[optional_loop_out],
2283 sum_outgoing += literal_lp_values[optional_loop_out];
2287 if (sum_outgoing < rhs_lower_bound - 1e-6) {
2288 manager->AddCut(outgoing.Build(),
"Circuit", lp_values);
2301 int num_nodes,
const std::vector<int>& tails,
const std::vector<int>& heads,
2302 const std::vector<Literal>& literals,
2306 if (num_nodes <= 2)
return;
2315 std::vector<Arc> relevant_arcs;
2318 std::vector<double> literal_lp_values(literals.size());
2319 std::vector<std::pair<double, int>> arc_by_decreasing_lp_values;
2321 for (
int i = 0; i < literals.size(); ++i) {
2323 const IntegerVariable direct_view = encoder->
GetLiteralView(literals[i]);
2325 lp_value = lp_values[direct_view];
2328 1.0 - lp_values[encoder->GetLiteralView(literals[i].Negated())];
2330 literal_lp_values[i] = lp_value;
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});
2336 std::sort(arc_by_decreasing_lp_values.begin(),
2337 arc_by_decreasing_lp_values.end(),
2338 std::greater<std::pair<double, int>>());
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) {
2355 auto get_root_and_compress_path = [&root](
int node) {
2357 while (root[r] != r) r = root[r];
2358 while (root[node] != r) {
2359 const int next = root[node];
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]);
2372 const int new_node = parent.size();
2373 parent.push_back(new_node);
2374 parent[
head] = new_node;
2375 parent[
tail] = new_node;
2379 root.push_back(new_node);
2380 root[
head] = new_node;
2381 root[
tail] = new_node;
2392 std::vector<int> pre_order(num_nodes);
2393 std::vector<absl::Span<const int>> subsets;
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);
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) {
2406 CHECK(graph[i].empty() || graph[i].size() == 2);
2407 if (parent[i] != i)
continue;
2412 while (!queue.empty()) {
2413 const int node = queue.back();
2418 const int start = start_index[node];
2419 if (new_size - start > 1) {
2420 subsets.emplace_back(&pre_order[start], new_size - start);
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);
2436 int64 total_demands = 0;
2437 if (!demands.empty()) {
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) {
2446 CHECK_LT(subset.size(), num_nodes);
2449 bool contain_depot =
false;
2450 int64 subset_demand = 0;
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];
2474 int64 min_outgoing_flow = 1;
2475 if (!demands.empty()) {
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;
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 min_outgoing_flow, lp_values, manager,
2513 for (
const int n : subset) in_subset[n] =
false;
2520 std::vector<IntegerVariable> GetAssociatedVariables(
2521 const std::vector<Literal>& literals, Model*
model) {
2522 auto* encoder =
model->GetOrCreate<IntegerEncoder>();
2523 std::vector<IntegerVariable> result;
2524 for (
const Literal l : literals) {
2525 const IntegerVariable direct_view = encoder->GetLiteralView(l);
2527 result.push_back(direct_view);
2529 result.push_back(encoder->GetLiteralView(l.Negated()));
2542 int num_nodes,
const std::vector<int>& tails,
const std::vector<int>& heads,
2543 const std::vector<Literal>& literals,
Model*
model) {
2545 result.
vars = GetAssociatedVariables(literals,
model);
2547 [num_nodes, tails, heads, literals,
model](
2551 num_nodes, tails, heads, literals, lp_values,
2552 {}, 0, manager,
model);
2558 const std::vector<int>& tails,
2559 const std::vector<int>& heads,
2560 const std::vector<Literal>& literals,
2561 const std::vector<int64>& demands,
2564 result.
vars = GetAssociatedVariables(literals,
model);
2566 [num_nodes, tails, heads, demands,
capacity, literals,
model](
2570 lp_values, demands,
capacity, manager,
2576 std::function<IntegerLiteral()>
2579 std::vector<IntegerVariable> variables;
2580 for (IntegerVariable
var : integer_variables_) {
2583 variables.push_back(
var);
2586 VLOG(1) <<
"HeuristicLPMostInfeasibleBinary has " << variables.size()
2589 return [
this, variables]() {
2593 double fractional_distance_best = -1.0;
2594 for (
const IntegerVariable
var : variables) {
2599 if (lb == ub)
continue;
2603 const double fractional_distance =
2605 lp_value - std::floor(lp_value +
kEpsilon));
2606 if (fractional_distance <
kEpsilon)
continue;
2609 if (fractional_distance > fractional_distance_best) {
2610 fractional_var =
var;
2611 fractional_distance_best = fractional_distance;
2625 std::vector<IntegerVariable> variables;
2626 for (IntegerVariable
var : integer_variables_) {
2629 variables.push_back(
var);
2632 VLOG(1) <<
"HeuristicLpReducedCostBinary has " << variables.size()
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);
2643 return [=]()
mutable {
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;
2657 for (
int i = 0; i < num_vars; i++) {
2658 const IntegerVariable
var = variables[i];
2663 if (lb == ub)
continue;
2667 if (std::abs(rc) <
kEpsilon)
continue;
2670 if (
value == 1.0 && rc < 0.0) {
2671 cost_to_zero[i] -= rc;
2672 num_cost_to_zero[i]++;
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];
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];
2692 if (selected_index >= 0) {
2700 void LinearProgrammingConstraint::UpdateAverageReducedCosts() {
2701 const int num_vars = integer_variables_.size();
2702 if (sum_cost_down_.size() < num_vars) {
2703 sum_cost_down_.resize(num_vars, 0.0);
2704 num_cost_down_.resize(num_vars, 0);
2705 sum_cost_up_.resize(num_vars, 0.0);
2706 num_cost_up_.resize(num_vars, 0);
2707 rc_scores_.resize(num_vars, 0.0);
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;
2719 num_calls_since_reduced_cost_averages_reset_ = 0;
2723 for (
int i = 0; i < num_vars; i++) {
2724 const IntegerVariable
var = integer_variables_[i];
2731 const double rc = lp_reduced_cost_[i];
2732 if (std::abs(rc) < kCpEpsilon)
continue;
2735 sum_cost_down_[i] -= rc;
2736 num_cost_down_[i]++;
2738 sum_cost_up_[i] += rc;
2745 rc_rev_int_repository_.
SetLevel(0);
2751 positions_by_decreasing_rc_score_.clear();
2752 for (
int i = 0; i < num_vars; i++) {
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);
2763 rc_scores_[i] = 0.5 * (a_down + a_up);
2768 if (rc_scores_[i] > 0.0) {
2769 positions_by_decreasing_rc_score_.push_back({-rc_scores_[i], i});
2772 std::sort(positions_by_decreasing_rc_score_.begin(),
2773 positions_by_decreasing_rc_score_.end());
2777 std::function<IntegerLiteral()>
2779 return [
this]() {
return this->LPReducedCostAverageDecision(); };
2782 IntegerLiteral LinearProgrammingConstraint::LPReducedCostAverageDecision() {
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];
2792 selected_index =
index;
2797 if (selected_index == -1)
return IntegerLiteral();
2798 const IntegerVariable
var = integer_variables_[selected_index];
2805 const IntegerValue value_ceil(
2807 if (value_ceil >= ub) {
2814 const IntegerValue value_floor(
2816 if (value_floor <= lb) {
2823 num_cost_up_[selected_index] > 0
2824 ? sum_cost_up_[selected_index] / num_cost_up_[selected_index]
2826 const double a_down =
2827 num_cost_down_[selected_index] > 0
2828 ? sum_cost_down_[selected_index] / num_cost_down_[selected_index]
2830 if (a_down < a_up) {
#define DCHECK_NE(val1, val2)
#define CHECK_LT(val1, val2)
#define CHECK_EQ(val1, val2)
#define CHECK_GE(val1, val2)
#define CHECK_GT(val1, val2)
#define CHECK_NE(val1, val2)
#define DCHECK_GT(val1, val2)
#define DCHECK(condition)
#define VLOG(verboselevel)
void assign(size_type n, const value_type &val)
void resize(size_type new_size)
void push_back(const value_type &x)
static int64 GCD64(int64 x, int64 y)
void SetLevel(int level) final
void SaveState(T *object)
A simple class to enforce both an elapsed time limit and a deterministic time limit in the same threa...
bool LimitReached()
Returns true when the external limit is true, or the deterministic time is over the deterministic lim...
void SetVariableBounds(ColIndex col, Fractional lower_bound, Fractional upper_bound)
void SetObjectiveOffset(Fractional objective_offset)
void SetCoefficient(RowIndex row, ColIndex col, Fractional value)
void SetConstraintBounds(RowIndex row, Fractional lower_bound, Fractional upper_bound)
ColIndex CreateNewVariable()
void AddSlackVariablesWhereNecessary(bool detect_integer_constraints)
void NotifyThatColumnsAreClean()
void SetObjectiveCoefficient(ColIndex col, Fractional value)
RowIndex CreateNewConstraint()
std::string GetDimensionString() const
Fractional objective_scaling_factor() const
const SparseColumn & GetSparseColumn(ColIndex col) const
RowIndex num_constraints() const
void Scale(LinearProgram *lp)
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)
const DenseRow & GetReducedCosts() const
int64 GetNumberOfIterations() const
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)
ProblemStatus GetProblemStatus() const
Fractional GetObjectiveValue() const
Fractional GetDualValue(RowIndex row) const
void ClearIntegralityScales()
void NotifyThatMatrixIsUnchangedForNextSolve()
ConstraintStatus GetConstraintStatus(RowIndex row) const
ColIndex GetProblemNumCols() const
void LoadStateForNextSolve(const BasisState &state)
RowIndex GetProblemNumRows() const
void ClearStateForNextSolve()
const BasisState & GetState() const
ColIndex GetBasis(RowIndex row) const
void SetParameters(const GlopParameters ¶meters)
const ScatteredRow & GetUnitRowLeftInverse(RowIndex row)
EntryIndex num_entries() const
LinearConstraint * mutable_cut()
bool TrySimpleKnapsack(const LinearConstraint base_ct, const std::vector< double > &lp_values, const std::vector< IntegerValue > &lower_bounds, const std::vector< IntegerValue > &upper_bounds)
void AlwaysCallAtLevelZero(int id)
void RegisterReversibleInt(int id, int *rev)
void WatchIntegerVariable(IntegerVariable i, int id, int watch_index=-1)
void WatchUpperBound(IntegerVariable var, int id, int watch_index=-1)
void SetPropagatorPriority(int id, int priority)
int Register(PropagatorInterface *propagator)
void AddLpVariable(IntegerVariable var)
void ProcessUpperBoundedConstraintWithSlackCreation(bool substitute_only_inner_variables, IntegerVariable first_slack, const absl::StrongVector< IntegerVariable, double > &lp_values, LinearConstraint *cut, std::vector< SlackInfo > *slack_infos)
bool DebugSlack(IntegerVariable first_slack, const LinearConstraint &initial_cut, const LinearConstraint &cut, const std::vector< SlackInfo > &info)
void SeparateSomeImpliedBoundCuts(const absl::StrongVector< IntegerVariable, double > &lp_values)
void AddData(double new_record)
double CurrentAverage() const
const IntegerVariable GetLiteralView(Literal lit) const
int NumLiftedBooleans() const
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)
ABSL_MUST_USE_RESULT bool Enqueue(IntegerLiteral i_lit, absl::Span< const Literal > literal_reason, absl::Span< const IntegerLiteral > integer_reason)
bool IsCurrentlyIgnored(IntegerVariable i) const
bool IsFixed(IntegerVariable i) const
IntegerLiteral LowerBoundAsLiteral(IntegerVariable i) const
bool ReportConflict(absl::Span< const Literal > literal_reason, absl::Span< const IntegerLiteral > integer_reason)
IntegerValue UpperBound(IntegerVariable i) const
IntegerValue LevelZeroUpperBound(IntegerVariable var) const
IntegerValue LevelZeroLowerBound(IntegerVariable var) const
void RelaxLinearReason(IntegerValue slack, absl::Span< const IntegerValue > coeffs, std::vector< IntegerLiteral > *reason) const
IntegerValue LowerBound(IntegerVariable i) const
IntegerLiteral UpperBoundAsLiteral(IntegerVariable i) const
bool IsFixedAtLevelZero(IntegerVariable var) const
void RemoveLevelZeroBounds(std::vector< IntegerLiteral > *reason) const
void RegisterReversibleClass(ReversibleInterface *rev)
bool ChangeLp(const absl::StrongVector< IntegerVariable, double > &lp_solution, glop::BasisState *solution_state)
bool DebugCheckConstraint(const LinearConstraint &cut)
void SetObjectiveCoefficient(IntegerVariable var, IntegerValue coeff)
ConstraintIndex Add(LinearConstraint ct, bool *added=nullptr)
const absl::StrongVector< ConstraintIndex, ConstraintInfo > & AllConstraints() const
const std::vector< ConstraintIndex > & LpConstraints() const
void AddAllConstraintsToLp()
bool AddCut(LinearConstraint ct, std::string type_name, const absl::StrongVector< IntegerVariable, double > &lp_solution, std::string extra_info="")
bool Propagate() override
double GetSolutionValue(IntegerVariable variable) const
void RegisterWith(Model *model)
glop::RowIndex ConstraintIndex
std::function< IntegerLiteral()> HeuristicLpReducedCostAverageBranching()
LinearProgrammingConstraint(Model *model)
std::function< IntegerLiteral()> HeuristicLpReducedCostBinary(Model *model)
~LinearProgrammingConstraint() override
void AddLinearConstraint(const LinearConstraint &ct)
bool IncrementalPropagate(const std::vector< int > &watch_indices) override
void SetLevel(int level) override
std::function< IntegerLiteral()> HeuristicLpMostInfeasibleBinary(Model *model)
void SetObjectiveCoefficient(IntegerVariable ivar, IntegerValue coeff)
double GetSolutionReducedCost(IntegerVariable variable) const
void AddCutGenerator(CutGenerator generator)
Class that owns everything related to a particular optimization model.
void ConvertToLinearConstraint(const std::vector< IntegerVariable > &integer_variables, IntegerValue upper_bound, LinearConstraint *result)
bool Add(glop::ColIndex col, IntegerValue value)
void ClearAndResize(int size)
std::vector< std::pair< glop::ColIndex, IntegerValue > > GetTerms()
bool AddLinearExpressionMultiple(IntegerValue multiplier, const std::vector< std::pair< glop::ColIndex, IntegerValue >> &terms)
void TransferToManager(const absl::StrongVector< IntegerVariable, double > &lp_solution, LinearConstraintManager *manager)
int CurrentDecisionLevel() const
void ProcessVariables(const std::vector< double > &lp_values, const std::vector< IntegerValue > &lower_bounds, const std::vector< IntegerValue > &upper_bounds)
std::vector< std::vector< std::pair< glop::RowIndex, IntegerValue > > > InterestingCandidates(ModelRandomGenerator *random)
void AddOneConstraint(glop::RowIndex, const std::vector< std::pair< glop::ColIndex, IntegerValue >> &terms, IntegerValue lb, IntegerValue ub)
static const int64 kint64max
static const int64 kint64min
const Collection::value_type::second_type & FindOrDie(const Collection &collection, const typename Collection::value_type::first_type &key)
StrictITIVector< ColIndex, Fractional > DenseRow
ColIndex RowToColIndex(RowIndex row)
RowIndex ColToRowIndex(ColIndex col)
StrictITIVector< RowIndex, Fractional > DenseColumn
IntegerValue FloorRatio(IntegerValue dividend, IntegerValue positive_divisor)
bool AddProductTo(IntegerValue a, IntegerValue b, IntegerValue *result)
constexpr IntegerValue kMaxIntegerValue(std::numeric_limits< IntegerValue::ValueType >::max() - 1)
IntType IntTypeAbs(IntType t)
void SeparateSubtourInequalities(int num_nodes, const std::vector< int > &tails, const std::vector< int > &heads, const std::vector< Literal > &literals, const absl::StrongVector< IntegerVariable, double > &lp_values, absl::Span< const int64 > demands, int64 capacity, LinearConstraintManager *manager, Model *model)
constexpr IntegerValue kMinIntegerValue(-kMaxIntegerValue)
const IntegerVariable kNoIntegerVariable(-1)
IntegerVariable PositiveVariable(IntegerVariable i)
CutGenerator CreateCVRPCutGenerator(int num_nodes, const std::vector< int > &tails, const std::vector< int > &heads, const std::vector< Literal > &literals, const std::vector< int64 > &demands, int64 capacity, Model *model)
std::vector< IntegerVariable > NegationOf(const std::vector< IntegerVariable > &vars)
IntegerValue ComputeInfinityNorm(const LinearConstraint &constraint)
bool VariableIsPositive(IntegerVariable i)
void DivideByGCD(LinearConstraint *constraint)
CutGenerator CreateStronglyConnectedGraphCutGenerator(int num_nodes, const std::vector< int > &tails, const std::vector< int > &heads, const std::vector< Literal > &literals, Model *model)
double ComputeActivity(const LinearConstraint &constraint, const absl::StrongVector< IntegerVariable, double > &values)
double ToDouble(IntegerValue value)
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
int64 CapSub(int64 x, int64 y)
std::pair< int64, int64 > Arc
int64 CapAdd(int64 x, int64 y)
int64 CapProd(int64 x, int64 y)
std::vector< IntegerVariable > vars
std::function< void(const absl::StrongVector< IntegerVariable, double > &lp_values, LinearConstraintManager *manager)> generate_cuts
static IntegerLiteral LowerOrEqual(IntegerVariable i, IntegerValue bound)
static IntegerLiteral GreaterOrEqual(IntegerVariable i, IntegerValue bound)
std::vector< IntegerValue > coeffs
std::vector< IntegerVariable > vars