24#include "absl/log/check.h"
25#include "absl/strings/str_cat.h"
26#include "absl/strings/string_view.h"
27#include "absl/types/span.h"
55using operations_research::MPConstraintProto;
56using operations_research::MPModelProto;
57using operations_research::MPVariableProto;
61void ScaleConstraint(absl::Span<const double> var_scaling,
64 for (
int i = 0;
i < num_terms; ++
i) {
65 const int var_index = mp_constraint->
var_index(
i);
71void ApplyVarScaling(absl::Span<const double> var_scaling,
74 for (
int i = 0;
i < num_variables; ++
i) {
75 const double scaling = var_scaling[
i];
87 ScaleConstraint(var_scaling, &mp_constraint);
89 for (MPGeneralConstraintProto& general_constraint :
91 switch (general_constraint.general_constraint_case()) {
93 ScaleConstraint(var_scaling,
94 general_constraint.mutable_indicator_constraint()
95 ->mutable_constraint());
103 LOG(FATAL) <<
"Scaling unsupported for general constraint of type "
104 << general_constraint.general_constraint_case();
114 std::vector<double> var_scaling(num_variables, 1.0);
115 for (
int i = 0;
i < num_variables; ++
i) {
117 if (max_bound == std::numeric_limits<double>::infinity()) {
118 var_scaling[
i] = scaling;
123 const double magnitude = std::max(std::abs(lb), std::abs(ub));
124 if (magnitude == 0 || magnitude > max_bound)
continue;
125 var_scaling[
i] = std::min(scaling, max_bound / magnitude);
127 ApplyVarScaling(var_scaling, mp_model);
135 const double initial_x = x;
138 int64_t current_q = 1;
140 while (current_q < limit) {
141 const double q =
static_cast<double>(current_q);
142 const double qx = q * initial_x;
143 const double qtolerance = q * tolerance;
144 if (std::abs(qx - std::round(qx)) < qtolerance) {
148 const double floored_x = std::floor(x);
149 if (floored_x >=
static_cast<double>(std::numeric_limits<int64_t>::max())) {
152 const int64_t new_q =
153 CapAdd(prev_q,
CapProd(
static_cast<int64_t
>(floored_x), current_q));
167double GetIntegralityMultiplier(
const MPModelProto& mp_model,
168 absl::Span<const double> var_scaling,
int var,
169 int ct_index,
double tolerance) {
172 double multiplier = 1.0;
173 double var_coeff = 0.0;
174 const double max_multiplier = 1e4;
188 if (multiplier == 0 || multiplier > max_multiplier)
return 0.0;
190 DCHECK_NE(var_coeff, 0.0);
194 if (!std::isfinite(bound))
continue;
195 if (std::abs(std::round(bound * multiplier) - bound * multiplier) >
196 tolerance * multiplier) {
200 return std::abs(multiplier * var_coeff);
210 int64_t num_changes = 0;
211 for (
int i = 0;
i < num_variables; ++
i) {
216 const double new_lb = std::isfinite(lb) ? std::ceil(lb - tolerance) : lb;
223 const double new_ub = std::isfinite(ub) ? std::floor(ub + tolerance) : ub;
229 if (new_ub < new_lb) {
230 SOLVER_LOG(logger,
"Empty domain for integer variable #",
i,
": [", lb,
241 int64_t num_variable_bounds_pushed_to_infinity = 0;
242 const double infinity = std::numeric_limits<double>::infinity();
243 for (
int i = 0;
i < num_variables; ++
i) {
246 if (std::isfinite(lb) && lb < -max_magnitude) {
247 ++num_variable_bounds_pushed_to_infinity;
252 if (std::isfinite(ub) && ub > max_magnitude) {
253 ++num_variable_bounds_pushed_to_infinity;
258 if (num_variable_bounds_pushed_to_infinity > 0) {
259 SOLVER_LOG(logger,
"Pushed ", num_variable_bounds_pushed_to_infinity,
260 " variable bounds to +/-infinity");
264 int64_t num_constraint_bounds_pushed_to_infinity = 0;
266 for (
int i = 0;
i < num_constraints; ++
i) {
269 if (std::isfinite(lb) && lb < -max_magnitude) {
270 ++num_constraint_bounds_pushed_to_infinity;
275 if (std::isfinite(ub) && ub > max_magnitude) {
276 ++num_constraint_bounds_pushed_to_infinity;
291 if (std::isfinite(lb) && lb < -max_magnitude) {
292 ++num_constraint_bounds_pushed_to_infinity;
297 if (std::isfinite(ub) && ub > max_magnitude) {
298 ++num_constraint_bounds_pushed_to_infinity;
303 if (num_constraint_bounds_pushed_to_infinity > 0) {
304 SOLVER_LOG(logger,
"Pushed ", num_constraint_bounds_pushed_to_infinity,
305 " constraint bounds to +/-infinity");
313 double max_dropped = 0.0;
316 for (
int i = 0;
i < num_variables; ++
i) {
320 max_dropped = std::max(max_dropped, std::abs(var->
lower_bound()));
325 max_dropped = std::max(max_dropped, std::abs(var->
upper_bound()));
330 for (
int i = 0;
i < num_constraints; ++
i) {
334 max_dropped = std::max(max_dropped, std::abs(ct->
lower_bound()));
339 max_dropped = std::max(max_dropped, std::abs(ct->
upper_bound()));
343 if (num_dropped > 0) {
344 SOLVER_LOG(logger,
"Set to zero ", num_dropped,
345 " variable or constraint bounds with largest magnitude ",
352 std::vector<double> max_bounds(num_variables);
353 for (
int i = 0;
i < num_variables; ++
i) {
357 max_bounds[
i] = value;
362 double largest_removed = 0.0;
369 int64_t num_removed = 0;
370 for (
int c = 0; c < num_constraints; ++c) {
374 if (size == 0)
continue;
375 const double threshold =
377 for (
int i = 0;
i < size; ++
i) {
380 if (std::abs(coeff) * max_bounds[var] < threshold) {
381 if (max_bounds[var] != 0) {
382 largest_removed = std::max(largest_removed, std::abs(coeff));
390 num_removed += size - new_size;
396 if (num_variables > 0) {
397 const double threshold =
399 for (
int var = 0; var < num_variables; ++var) {
401 if (coeff == 0.0)
continue;
402 if (std::abs(coeff) * max_bounds[var] < threshold) {
404 if (max_bounds[var] != 0) {
405 largest_removed = std::max(largest_removed, std::abs(coeff));
412 if (num_removed > 0) {
414 " near zero terms with largest magnitude of ", largest_removed,
425 switch (general_constraint.general_constraint_case()) {
433 SOLVER_LOG(logger,
"General constraints of type ",
434 general_constraint.general_constraint_case(),
435 " are not supported.");
443 for (
int i = 0;
i < num_variables; ++
i) {
454 SOLVER_LOG(logger,
"Objective coefficient is too large: ",
462 for (
int c = 0; c < num_constraints; ++c) {
473 if (std::abs(coeff) > threshold) {
474 SOLVER_LOG(logger,
"Constraint coefficient is too large: ", coeff);
486 std::vector<double> var_scaling(num_variables, 1.0);
488 int initial_num_integers = 0;
489 for (
int i = 0;
i < num_variables; ++
i) {
492 VLOG(1) <<
"Initial num integers: " << initial_num_integers;
495 const double tolerance = 1e-6;
496 std::vector<int> constraint_queue;
499 std::vector<int> constraint_to_num_non_integer(num_constraints, 0);
500 std::vector<std::vector<int>> var_to_constraints(num_variables);
501 for (
int i = 0;
i < num_constraints; ++
i) {
504 for (
const int var : mp_constraint.
var_index()) {
506 var_to_constraints[var].push_back(
i);
507 constraint_to_num_non_integer[
i]++;
510 if (constraint_to_num_non_integer[
i] == 1) {
511 constraint_queue.push_back(
i);
514 VLOG(1) <<
"Initial constraint queue: " << constraint_queue.size() <<
" / "
517 int num_detected = 0;
518 double max_scaling = 0.0;
519 auto scale_and_mark_as_integer = [&](
int var,
double scaling)
mutable {
522 CHECK_EQ(var_scaling[var], 1.0);
523 if (scaling != 1.0) {
524 VLOG(2) <<
"Scaled " << var <<
" by " << scaling;
528 max_scaling = std::max(max_scaling, scaling);
532 var_scaling[var] = scaling;
536 for (
const int ct_index : var_to_constraints[var]) {
537 constraint_to_num_non_integer[ct_index]--;
538 if (constraint_to_num_non_integer[ct_index] == 1) {
539 constraint_queue.push_back(ct_index);
544 int num_fail_due_to_rhs = 0;
545 int num_fail_due_to_large_multiplier = 0;
546 int num_processed_constraints = 0;
547 while (!constraint_queue.empty()) {
548 const int top_ct_index = constraint_queue.back();
549 constraint_queue.pop_back();
553 if (constraint_to_num_non_integer[top_ct_index] == 0)
continue;
559 ++num_processed_constraints;
571 double multiplier = 1.0;
572 const double max_multiplier = 1e4;
586 if (multiplier == 0 || multiplier > max_multiplier) {
592 if (multiplier == 0 || multiplier > max_multiplier) {
593 ++num_fail_due_to_large_multiplier;
599 if (std::abs(std::round(rhs * multiplier) - rhs * multiplier) >
600 tolerance * multiplier) {
601 ++num_fail_due_to_rhs;
611 double best_scaling = std::abs(var_coeff * multiplier);
612 for (
const int ct_index : var_to_constraints[var]) {
613 if (ct_index == top_ct_index)
continue;
614 if (constraint_to_num_non_integer[ct_index] != 1)
continue;
620 const double multiplier = GetIntegralityMultiplier(
621 *mp_model, var_scaling, var, ct_index, tolerance);
622 if (multiplier != 0.0 && multiplier < best_scaling) {
623 best_scaling = multiplier;
627 scale_and_mark_as_integer(var, best_scaling);
635 int num_in_inequalities = 0;
636 int num_to_be_handled = 0;
637 for (
int var = 0; var < num_variables; ++var) {
641 if (var_to_constraints[var].empty())
continue;
644 for (
const int ct_index : var_to_constraints[var]) {
645 if (constraint_to_num_non_integer[ct_index] != 1) {
652 std::vector<double> scaled_coeffs;
653 for (
const int ct_index : var_to_constraints[var]) {
654 const double multiplier = GetIntegralityMultiplier(
655 *mp_model, var_scaling, var, ct_index, tolerance);
656 if (multiplier == 0.0) {
660 scaled_coeffs.push_back(multiplier);
669 double scaling = scaled_coeffs[0];
670 for (
const double c : scaled_coeffs) {
671 scaling = std::min(scaling, c);
673 CHECK_GT(scaling, 0.0);
674 for (
const double c : scaled_coeffs) {
675 const double fraction = c / scaling;
676 if (std::abs(std::round(fraction) - fraction) > tolerance) {
690 if (!std::isfinite(bound))
continue;
691 if (std::abs(std::round(bound * scaling) - bound * scaling) >
692 tolerance * scaling) {
704 ++num_in_inequalities;
705 scale_and_mark_as_integer(var, scaling);
707 VLOG(1) <<
"num_new_integer: " << num_detected
708 <<
" num_processed_constraints: " << num_processed_constraints
709 <<
" num_rhs_fail: " << num_fail_due_to_rhs
710 <<
" num_multiplier_fail: " << num_fail_due_to_large_multiplier;
712 if (num_to_be_handled > 0) {
713 SOLVER_LOG(logger,
"Missed ", num_to_be_handled,
714 " potential implied integer.");
717 const int num_integers = initial_num_integers + num_detected;
718 SOLVER_LOG(logger,
"Num integers: ", num_integers,
"/", num_variables,
719 " (implied: ", num_detected,
720 " in_inequalities: ", num_in_inequalities,
721 " max_scaling: ", max_scaling,
")",
722 (num_integers == num_variables ?
" [IP] " :
" [MIP] "));
724 ApplyVarScaling(var_scaling, mp_model);
731struct ConstraintScaler {
733 ConstraintProto* AddConstraint(
const MPModelProto& mp_model,
735 CpModelProto* cp_model);
737 bool keep_names =
false;
738 double max_relative_coeff_error = 0.0;
739 double max_absolute_rhs_error = 0.0;
740 double max_scaling_factor = 0.0;
741 double min_scaling_factor = std::numeric_limits<double>::infinity();
743 double wanted_precision = 1e-6;
744 int64_t scaling_target = int64_t{1} << 50;
745 std::vector<int> var_indices;
746 std::vector<double> coefficients;
747 std::vector<double> lower_bounds;
748 std::vector<double> upper_bounds;
759 auto* constraint = cp_model->add_constraints();
760 if (keep_names) constraint->
set_name(mp_constraint.
name());
761 auto* arg = constraint->mutable_linear();
766 coefficients.clear();
767 lower_bounds.clear();
768 upper_bounds.clear();
770 for (
int i = 0;
i < num_coeffs; ++
i) {
771 const auto& var_proto = cp_model->variables(mp_constraint.
var_index(i));
772 const int64_t lb = var_proto.domain(0);
773 const int64_t ub = var_proto.domain(var_proto.domain_size() - 1);
774 if (lb == 0 && ub == 0)
continue;
777 if (coeff == 0.0)
continue;
779 var_indices.push_back(mp_constraint.
var_index(i));
780 coefficients.push_back(coeff);
781 lower_bounds.push_back(lb);
782 upper_bounds.push_back(ub);
785 double relative_coeff_error;
786 double scaled_sum_error;
788 coefficients, lower_bounds, upper_bounds, scaling_target,
789 wanted_precision, &relative_coeff_error, &scaled_sum_error);
790 if (scaling_factor == 0.0) {
794 LOG(DFATAL) <<
"Scaling factor of zero while scaling constraint: "
800 max_relative_coeff_error =
801 std::max(relative_coeff_error, max_relative_coeff_error);
802 max_scaling_factor = std::max(scaling_factor / gcd, max_scaling_factor);
803 min_scaling_factor = std::min(scaling_factor / gcd, min_scaling_factor);
805 for (
int i = 0;
i < coefficients.size(); ++
i) {
806 const double scaled_value = coefficients[
i] * scaling_factor;
807 const int64_t value =
static_cast<int64_t
>(std::round(scaled_value)) / gcd;
809 arg->add_vars(var_indices[i]);
810 arg->add_coeffs(value);
813 max_absolute_rhs_error =
814 std::max(max_absolute_rhs_error, scaled_sum_error / scaling_factor);
830 const Fractional scaled_lb = std::ceil(lb * scaling_factor);
831 if (lb == kInfinity || scaled_lb >= std::numeric_limits<int64_t>::max()) {
833 arg->add_domain(std::numeric_limits<int64_t>::max());
834 }
else if (lb == -kInfinity ||
835 scaled_lb <= std::numeric_limits<int64_t>::min()) {
836 arg->add_domain(std::numeric_limits<int64_t>::min());
838 arg->add_domain(
CeilRatio(IntegerValue(
static_cast<int64_t
>(scaled_lb)),
843 const Fractional scaled_ub = std::floor(ub * scaling_factor);
844 if (ub == -kInfinity || scaled_ub <= std::numeric_limits<int64_t>::min()) {
846 arg->add_domain(std::numeric_limits<int64_t>::min());
847 }
else if (ub == kInfinity ||
848 scaled_ub >= std::numeric_limits<int64_t>::max()) {
849 arg->add_domain(std::numeric_limits<int64_t>::max());
851 arg->add_domain(
FloorRatio(IntegerValue(
static_cast<int64_t
>(scaled_ub)),
860double FindFractionalScaling(absl::Span<const double> coefficients,
862 double multiplier = 1.0;
863 for (
const double coeff : coefficients) {
865 multiplier * tolerance);
866 if (multiplier == 0.0)
break;
874 absl::Span<const double> coefficients,
875 absl::Span<const double> lower_bounds,
876 absl::Span<const double> upper_bounds, int64_t max_absolute_activity,
877 double wanted_absolute_activity_precision,
double* relative_coeff_error,
878 double* scaled_sum_error) {
881 coefficients, lower_bounds, upper_bounds, max_absolute_activity);
882 if (scaling_factor == 0.0)
return scaling_factor;
892 double x = std::min(scaling_factor, 1.0);
893 for (; x <= scaling_factor; x *= 2) {
895 relative_coeff_error, scaled_sum_error);
896 if (*scaled_sum_error < wanted_absolute_activity_precision * x)
break;
899 if (x == scaling_factor)
break;
902 DCHECK(std::isfinite(scaling_factor));
912 const double integer_factor = FindFractionalScaling(coefficients, 1e-8);
913 DCHECK(std::isfinite(integer_factor));
914 if (integer_factor != 0 && integer_factor < scaling_factor) {
915 double local_relative_coeff_error;
916 double local_scaled_sum_error;
918 integer_factor, &local_relative_coeff_error,
919 &local_scaled_sum_error);
920 if (local_scaled_sum_error * scaling_factor <=
921 *scaled_sum_error * integer_factor ||
922 local_scaled_sum_error <
923 wanted_absolute_activity_precision * integer_factor) {
924 *relative_coeff_error = local_relative_coeff_error;
925 *scaled_sum_error = local_scaled_sum_error;
926 scaling_factor = integer_factor;
930 DCHECK(std::isfinite(scaling_factor));
931 return scaling_factor;
938 CHECK(cp_model !=
nullptr);
954 const int64_t kMaxVariableBound =
957 int num_truncated_bounds = 0;
958 int num_small_domains = 0;
959 const int64_t kSmallDomainSize = 1000;
965 for (
int i = 0;
i < num_variables; ++
i) {
976 if (mp_var.
lower_bound() >
static_cast<double>(kMaxVariableBound) ||
977 mp_var.
upper_bound() <
static_cast<double>(-kMaxVariableBound)) {
978 SOLVER_LOG(logger,
"Error: variable ", mp_var,
979 " is outside [-mip_max_bound..mip_max_bound]");
984 for (
const bool lower : {
true,
false}) {
986 if (std::abs(bound) + kWantedPrecision >=
987 static_cast<double>(kMaxVariableBound)) {
988 ++num_truncated_bounds;
989 cp_var->
add_domain(bound < 0 ? -kMaxVariableBound : kMaxVariableBound);
995 static_cast<int64_t
>(lower ? std::ceil(bound - kWantedPrecision)
996 : std::floor(bound + kWantedPrecision)));
1000 LOG(WARNING) <<
"Variable #" <<
i <<
" cannot take integer value. "
1009 if (diff > kWantedPrecision && diff < kSmallDomainSize) {
1010 ++num_small_domains;
1015 if (num_truncated_bounds > 0) {
1016 SOLVER_LOG(logger,
"Warning: ", num_truncated_bounds,
1017 " bounds were truncated to ", kMaxVariableBound,
".");
1019 if (num_small_domains > 0) {
1020 SOLVER_LOG(logger,
"Warning: ", num_small_domains,
1021 " continuous variable domain with fewer than ", kSmallDomainSize,
1025 ConstraintScaler scaler;
1026 const int64_t kScalingTarget = int64_t{1}
1028 scaler.wanted_precision = kWantedPrecision;
1029 scaler.scaling_target = kScalingTarget;
1030 scaler.keep_names = keep_names;
1034 scaler.AddConstraint(mp_model, mp_constraint, cp_model);
1038 switch (general_constraint.general_constraint_case()) {
1040 const auto& indicator_constraint =
1041 general_constraint.indicator_constraint();
1043 indicator_constraint.constraint();
1045 scaler.AddConstraint(mp_model, mp_constraint, cp_model);
1046 if (ct ==
nullptr)
continue;
1049 const int var = indicator_constraint.var_index();
1050 const int value = indicator_constraint.var_value();
1055 const auto& and_constraint = general_constraint.and_constraint();
1056 absl::string_view name = general_constraint.name();
1059 ct_pos->
set_name(name.empty() ?
"" : absl::StrCat(name,
"_pos"));
1062 and_constraint.var_index();
1065 ct_neg->
set_name(name.empty() ?
"" : absl::StrCat(name,
"_neg"));
1067 NegatedRef(and_constraint.resultant_var_index()));
1068 for (
const int var_index : and_constraint.var_index()) {
1074 const auto& or_constraint = general_constraint.or_constraint();
1075 absl::string_view name = general_constraint.name();
1078 ct_pos->
set_name(name.empty() ?
"" : absl::StrCat(name,
"_pos"));
1081 or_constraint.var_index();
1084 ct_neg->
set_name(name.empty() ?
"" : absl::StrCat(name,
"_neg"));
1086 NegatedRef(or_constraint.resultant_var_index()));
1087 for (
const int var_index : or_constraint.var_index()) {
1093 LOG(ERROR) <<
"Can't convert general constraints of type "
1094 << general_constraint.general_constraint_case()
1095 <<
" to CpModelProto.";
1101 SOLVER_LOG(logger,
"Maximum constraint coefficient relative error: ",
1102 scaler.max_relative_coeff_error);
1103 SOLVER_LOG(logger,
"Maximum constraint worst-case activity error: ",
1104 scaler.max_absolute_rhs_error,
1106 ?
" [Potentially IMPRECISE]"
1108 SOLVER_LOG(logger,
"Constraint scaling factor range: [",
1109 scaler.min_scaling_factor,
", ", scaler.max_scaling_factor,
"]");
1117 for (
int i = 0;
i < num_variables; ++
i) {
1120 float_objective->add_vars(
i);
1126 if (float_objective->offset() == 0 && float_objective->vars().empty()) {
1136 for (
const int ref : literals) {
1153 CHECK(output !=
nullptr);
1157 const int num_vars =
input.variables().size();
1158 for (
int v = 0; v < num_vars; ++v) {
1159 if (
input.variables(v).domain().size() != 2) {
1160 VLOG(1) <<
"Cannot convert "
1172 if (
input.has_objective()) {
1173 double factor =
input.objective().scaling_factor();
1174 if (factor == 0.0) factor = 1.0;
1175 const int num_terms =
input.objective().vars().size();
1176 for (
int i = 0;
i < num_terms; ++
i) {
1177 const int var =
input.objective().vars(
i);
1178 if (var < 0)
return false;
1181 factor *
input.objective().coeffs(
i));
1187 }
else if (
input.has_floating_point_objective()) {
1188 const int num_terms =
input.floating_point_objective().vars().size();
1189 for (
int i = 0;
i < num_terms; ++
i) {
1190 const int var =
input.floating_point_objective().vars(
i);
1191 if (var < 0)
return false;
1194 input.floating_point_objective().coeffs(
i));
1203 const int num_constraints =
input.constraints().size();
1204 std::vector<int> tmp_literals;
1205 for (
int c = 0; c < num_constraints; ++c) {
1237 tmp_literals.clear();
1243 tmp_literals.push_back(ref);
1244 const int shift = AppendSumOfLiteral(tmp_literals, out);
1247 tmp_literals.pop_back();
1253 VLOG(1) <<
"Cannot convert constraint: "
1259 int64_t min_activity = 0;
1260 int64_t max_activity = 0;
1261 const int num_terms = ct.
linear().
vars().size();
1262 for (
int i = 0;
i < num_terms; ++
i) {
1264 if (var < 0)
return false;
1265 DCHECK_EQ(
input.variables(var).domain().size(), 2);
1268 min_activity += coeff *
input.variables(var).domain(0);
1269 max_activity += coeff *
input.variables(var).domain(1);
1271 min_activity += coeff *
input.variables(var).domain(1);
1272 max_activity += coeff *
input.variables(var).domain(0);
1288 for (
int i = 0;
i < num_terms; ++
i) {
1290 if (var < 0)
return false;
1297 std::vector<MPConstraintProto*> out_cts;
1302 const int64_t coeff = max_activity - ct.
linear().
domain(1);
1315 out_cts.push_back(high_out_ct);
1334 out_cts.push_back(low_out_ct);
1337 for (
int i = 0;
i < num_terms; ++
i) {
1339 if (var < 0)
return false;
1340 out_ct->add_var_index(var);
1356 absl::Span<
const std::pair<int, double>> objective,
1357 double objective_offset,
bool maximize,
1363 std::vector<int> var_indices;
1364 std::vector<double> coefficients;
1365 std::vector<double> lower_bounds;
1366 std::vector<double> upper_bounds;
1367 double min_magnitude = std::numeric_limits<double>::infinity();
1368 double max_magnitude = 0.0;
1369 double l1_norm = 0.0;
1370 for (
const auto& [var, coeff] : objective) {
1371 const auto& var_proto = cp_model->
variables(var);
1372 const int64_t lb = var_proto.
domain(0);
1373 const int64_t ub = var_proto.domain(var_proto.domain_size() - 1);
1375 if (lb != 0) objective_offset += lb * coeff;
1378 var_indices.push_back(var);
1379 coefficients.push_back(coeff);
1380 lower_bounds.push_back(lb);
1381 upper_bounds.push_back(ub);
1383 min_magnitude = std::min(min_magnitude, std::abs(coeff));
1384 max_magnitude = std::max(max_magnitude, std::abs(coeff));
1385 l1_norm += std::abs(coeff);
1388 if (coefficients.empty() && objective_offset == 0.0)
return true;
1390 if (!coefficients.empty()) {
1391 const double average_magnitude =
1392 l1_norm /
static_cast<double>(coefficients.size());
1393 SOLVER_LOG(logger,
"[Scaling] Floating point objective has ",
1394 coefficients.size(),
" terms with magnitude in [", min_magnitude,
1395 ", ", max_magnitude,
"] average = ", average_magnitude);
1399 const int64_t max_absolute_activity = int64_t{1}
1401 const double wanted_precision =
1404 double relative_coeff_error;
1405 double scaled_sum_error;
1407 coefficients, lower_bounds, upper_bounds, max_absolute_activity,
1408 wanted_precision, &relative_coeff_error, &scaled_sum_error);
1409 if (scaling_factor == 0.0) {
1410 LOG(ERROR) <<
"Scaling factor of zero while scaling objective! This "
1411 "likely indicate an infinite coefficient in the objective.";
1418 SOLVER_LOG(logger,
"[Scaling] Objective coefficient relative error: ",
1419 relative_coeff_error);
1420 SOLVER_LOG(logger,
"[Scaling] Objective worst-case absolute error: ",
1421 scaled_sum_error / scaling_factor);
1423 "[Scaling] Objective scaling factor: ", scaling_factor / gcd);
1425 if (scaled_sum_error / scaling_factor > wanted_precision) {
1427 "[Scaling] Warning: the worst-case absolute error is greater "
1428 "than the wanted precision (",
1430 "). Try to increase mip_max_activity_exponent (default = ",
1432 ") or reduced your variables range and/or objective "
1433 "coefficient. We will continue the solve, but the final "
1434 "objective value might be off.");
1441 const int64_t mult = maximize ? -1 : 1;
1442 objective_proto->
set_offset(objective_offset * scaling_factor / gcd * mult);
1443 objective_proto->set_scaling_factor(1.0 / scaling_factor * gcd * mult);
1444 for (
int i = 0;
i < coefficients.size(); ++
i) {
1445 const int64_t value =
1446 static_cast<int64_t
>(std::round(coefficients[
i] * scaling_factor)) /
1449 objective_proto->add_vars(var_indices[
i]);
1450 objective_proto->add_coeffs(value * mult);
1454 if (scaled_sum_error == 0.0) {
1455 objective_proto->set_scaling_was_exact(
true);
1463 CHECK(problem !=
nullptr);
1471 for (
int var_id(0); var_id < num_variables; ++var_id) {
1482 if (lb <= -1.0) is_binary =
false;
1483 if (ub >= 2.0) is_binary =
false;
1486 if (lb <= 0.0 && ub >= 1.0) {
1488 }
else if (lb <= 1.0 && ub >= 1.0) {
1495 }
else if (lb <= 0.0 && ub >= 0.0) {
1510 LOG(WARNING) <<
"The variable #" << var_id <<
" with name "
1511 << mp_var.
name() <<
" is not binary. "
1512 <<
"lb: " << lb <<
" ub: " << ub;
1518 const int64_t kInt64Max = std::numeric_limits<int64_t>::max();
1519 double max_relative_error = 0.0;
1520 double max_bound_error = 0.0;
1521 double max_scaling_factor = 0.0;
1522 double relative_error = 0.0;
1523 double scaling_factor = 0.0;
1524 std::vector<double> coefficients;
1532 coefficients.clear();
1534 for (
int i = 0;
i < num_coeffs; ++
i) {
1541 max_relative_error = std::max(relative_error, max_relative_error);
1542 max_scaling_factor = std::max(scaling_factor / gcd, max_scaling_factor);
1544 double bound_error = 0.0;
1545 for (
int i = 0;
i < num_coeffs; ++
i) {
1546 const double scaled_value = mp_constraint.
coefficient(
i) * scaling_factor;
1547 bound_error += std::abs(round(scaled_value) - scaled_value);
1548 const int64_t value =
static_cast<int64_t
>(round(scaled_value)) / gcd;
1554 max_bound_error = std::max(max_bound_error, bound_error);
1563 if (lb * scaling_factor >
static_cast<double>(kInt64Max)) {
1564 LOG(WARNING) <<
"A constraint is trivially unsatisfiable.";
1567 if (lb * scaling_factor > -
static_cast<double>(kInt64Max)) {
1570 static_cast<int64_t
>(round(lb * scaling_factor - bound_error)) /
1576 if (ub * scaling_factor < -
static_cast<double>(kInt64Max)) {
1577 LOG(WARNING) <<
"A constraint is trivially unsatisfiable.";
1580 if (ub * scaling_factor <
static_cast<double>(kInt64Max)) {
1583 static_cast<int64_t
>(round(ub * scaling_factor + bound_error)) /
1590 LOG(INFO) <<
"Maximum constraint relative error: " << max_relative_error;
1591 LOG(INFO) <<
"Maximum constraint bound error: " << max_bound_error;
1592 LOG(INFO) <<
"Maximum constraint scaling factor: " << max_scaling_factor;
1595 coefficients.clear();
1596 for (
int var_id = 0; var_id < num_variables; ++var_id) {
1603 max_relative_error = std::max(relative_error, max_relative_error);
1606 LOG(INFO) <<
"objective relative error: " << relative_error;
1607 LOG(INFO) <<
"objective scaling factor: " << scaling_factor / gcd;
1615 for (
int var_id = 0; var_id < num_variables; ++var_id) {
1617 const int64_t value =
1618 static_cast<int64_t
>(
1631 const double kRelativeTolerance = 1e-8;
1632 if (max_relative_error > kRelativeTolerance) {
1633 LOG(WARNING) <<
"The relative error during double -> int64_t conversion "
1661 for (
int i = 0;
i < constraint.literals_size(); ++
i) {
1662 const int literal = constraint.literals(
i);
1663 const double coeff = constraint.coefficients(
i);
1664 const ColIndex variable_index = ColIndex(abs(literal) - 1);
1674 constraint.has_lower_bound() ? constraint.lower_bound() - sum
1676 constraint.has_upper_bound() ? constraint.upper_bound() - sum
1686 const int literal = objective.
literals(
i);
1687 const double coeff =
1688 static_cast<double>(objective.
coefficients(
i)) * scaling_factor;
1689 const ColIndex variable_index = ColIndex(abs(literal) - 1);
1705 const CpModelProto& model_proto_with_floating_point_objective,
1707 const int64_t inner_integer_objective_lower_bound) {
1710 const CpModelProto& proto = model_proto_with_floating_point_objective;
1714 static_cast<double>(domain[domain.size() - 1]));
1722 for (
int i = 0;
i < float_obj.
vars().size(); ++
i) {
1723 const glop::ColIndex col(float_obj.
vars(
i));
1731 ct,
static_cast<double>(inner_integer_objective_lower_bound),
1732 std::numeric_limits<double>::infinity());
1733 for (
int i = 0;
i < integer_objective.
vars().size(); ++
i) {
1735 static_cast<double>(integer_objective.
coeffs(
i)));
1753 return float_obj.
maximize() ? std::numeric_limits<double>::infinity()
1754 : -std::numeric_limits<double>::infinity();
void set_var_index(int index, ::int32_t value)
::google::protobuf::RepeatedField< double > *PROTOBUF_NONNULL mutable_coefficient()
double lower_bound() const
::google::protobuf::RepeatedField<::int32_t > *PROTOBUF_NONNULL mutable_var_index()
void set_lower_bound(double value)
void add_coefficient(double value)
void set_upper_bound(double value)
double coefficient(int index) const
void set_coefficient(int index, double value)
void add_var_index(::int32_t value)
::int32_t var_index(int index) const
double upper_bound() const
GeneralConstraintCase general_constraint_case() const
::operations_research::MPIndicatorConstraint *PROTOBUF_NONNULL mutable_indicator_constraint()
::operations_research::MPConstraintProto *PROTOBUF_NONNULL mutable_constraint()
const ::operations_research::MPVariableProto & variable(int index) const
int variable_size() const
repeated .operations_research.MPVariableProto variable = 3;
void set_maximize(bool value)
const ::operations_research::MPGeneralConstraintProto & general_constraint(int index) const
::operations_research::MPConstraintProto *PROTOBUF_NONNULL add_constraint()
void clear_objective_offset()
::operations_research::MPConstraintProto *PROTOBUF_NONNULL mutable_constraint(int index)
::operations_research::MPVariableProto *PROTOBUF_NONNULL add_variable()
const ::operations_research::MPConstraintProto & constraint(int index) const
void set_objective_offset(double value)
int general_constraint_size() const
repeated .operations_research.MPGeneralConstraintProto general_constraint = 7;
ABSL_ATTRIBUTE_REINITIALIZES void Clear() PROTOBUF_FINAL
::operations_research::MPVariableProto *PROTOBUF_NONNULL mutable_variable(int index)
const ::std::string & name() const
::operations_research::MPGeneralConstraintProto *PROTOBUF_NONNULL mutable_general_constraint(int index)
int constraint_size() const
repeated .operations_research.MPConstraintProto constraint = 4;
double objective_offset() const
void set_upper_bound(double value)
void set_is_integer(bool value)
void set_objective_coefficient(double value)
double lower_bound() const
double upper_bound() const
void clear_objective_coefficient()
void set_lower_bound(double value)
const ::std::string & name() const
double objective_coefficient() const
void set_change_status_to_imprecise(bool value)
void set_max_number_of_iterations(::int64_t value)
A full-fledged linear programming solver.
Fractional GetObjectiveValue() const
Returns the objective value of the solution with its offset and scaling.
void SetParameters(const GlopParameters ¶meters)
ABSL_MUST_USE_RESULT ProblemStatus Solve(const LinearProgram &lp)
const DenseRow & objective_coefficients() const
Returns the objective coefficients (or cost) of variables as a row vector.
void Clear()
Clears, i.e. reset the object to its initial value.
void SetConstraintName(RowIndex row, absl::string_view name)
@ INTEGER
The variable must only take integer values.
void SetObjectiveOffset(Fractional objective_offset)
void SetObjectiveCoefficient(ColIndex col, Fractional value)
void SetVariableBounds(ColIndex col, Fractional lower_bound, Fractional upper_bound)
void SetVariableType(ColIndex col, VariableType type)
Set the type of the variable.
void SetConstraintBounds(RowIndex row, Fractional lower_bound, Fractional upper_bound)
void SetCoefficient(RowIndex row, ColIndex col, Fractional value)
Defines the coefficient for col / row.
void SetVariableName(ColIndex col, absl::string_view name)
RowIndex CreateNewConstraint()
ColIndex CreateNewVariable()
void SetMaximizationProblem(bool maximize)
::int32_t literals(int index) const
void add_literals(::int32_t value)
::google::protobuf::RepeatedField<::int32_t > *PROTOBUF_NONNULL mutable_literals()
ConstraintCase constraint_case() const
::operations_research::sat::BoolArgumentProto *PROTOBUF_NONNULL mutable_bool_and()
::int32_t enforcement_literal(int index) const
const ::operations_research::sat::BoolArgumentProto & exactly_one() const
const ::operations_research::sat::LinearConstraintProto & linear() const
const ::operations_research::sat::BoolArgumentProto & at_most_one() const
const ::operations_research::sat::BoolArgumentProto & bool_and() const
::operations_research::sat::BoolArgumentProto *PROTOBUF_NONNULL mutable_bool_or()
void add_enforcement_literal(::int32_t value)
const ::operations_research::sat::BoolArgumentProto & bool_or() const
void set_name(Arg_ &&arg, Args_... args)
void clear_floating_point_objective()
::operations_research::sat::FloatObjectiveProto *PROTOBUF_NONNULL mutable_floating_point_objective()
const ::operations_research::sat::IntegerVariableProto & variables(int index) const
const ::operations_research::sat::FloatObjectiveProto & floating_point_objective() const
ABSL_ATTRIBUTE_REINITIALIZES void Clear() PROTOBUF_FINAL
void set_name(Arg_ &&arg, Args_... args)
::operations_research::sat::IntegerVariableProto *PROTOBUF_NONNULL add_variables()
::operations_research::sat::ConstraintProto *PROTOBUF_NONNULL add_constraints()
::operations_research::sat::CpObjectiveProto *PROTOBUF_NONNULL mutable_objective()
void set_offset(double value)
::int32_t vars(int index) const
::int64_t coeffs(int index) const
::int32_t vars(int index) const
void set_maximize(bool value)
double coeffs(int index) const
void add_domain(::int64_t value)
void set_name(Arg_ &&arg, Args_... args)
::int64_t domain(int index) const
void set_upper_bound(::int64_t value)
void set_name(Arg_ &&arg, Args_... args)
void set_lower_bound(::int64_t value)
void add_literals(::int32_t value)
void add_coefficients(::int64_t value)
void set_name(Arg_ &&arg, Args_... args)
::operations_research::sat::LinearBooleanConstraint *PROTOBUF_NONNULL add_constraints()
const ::operations_research::sat::LinearBooleanConstraint & constraints(int index) const
const ::operations_research::sat::LinearObjective & objective() const
void set_num_variables(::int32_t value)
::std::string *PROTOBUF_NONNULL add_var_names()
const ::std::string & var_names(int index) const
ABSL_ATTRIBUTE_REINITIALIZES void Clear() PROTOBUF_FINAL
::int32_t num_variables() const
int var_names_size() const
repeated string var_names = 6;
::operations_research::sat::LinearObjective *PROTOBUF_NONNULL mutable_objective()
::int32_t vars(int index) const
::int64_t coeffs(int index) const
::int64_t domain(int index) const
int literals_size() const
repeated int32 literals = 1;
void add_literals(::int32_t value)
void add_coefficients(::int64_t value)
::int64_t coefficients(int index) const
void set_scaling_factor(double value)
double scaling_factor() const
::int32_t literals(int index) const
void set_offset(double value)
double lower_bound() const
const ::std::string & name() const
int coefficient_size() const
repeated double coefficient = 7 [packed = true];
double coefficient(int index) const
void set_coefficient(int index, double value)
::int32_t var_index(int index) const
double upper_bound() const
const ::operations_research::MPVariableProto & variable(int index) const
int variable_size() const
repeated .operations_research.MPVariableProto variable = 3;
::operations_research::MPConstraintProto *PROTOBUF_NONNULL mutable_constraint(int index)
::operations_research::MPVariableProto *PROTOBUF_NONNULL mutable_variable(int index)
::operations_research::MPGeneralConstraintProto *PROTOBUF_NONNULL mutable_general_constraint(int index)
double lower_bound() const
double upper_bound() const
double objective_coefficient() const
double mip_check_precision() const
double mip_max_valid_magnitude() const
bool ignore_names() const
::int32_t mip_max_activity_exponent() const
double mip_wanted_precision() const
double mip_max_bound() const
double mip_drop_tolerance() const
double absolute_gap_limit() const
constexpr Fractional kInfinity
Infinity for type Fractional.
ProblemStatus
Different statuses for a given problem.
IntegerValue FloorRatio(IntegerValue dividend, IntegerValue positive_divisor)
double FindBestScalingAndComputeErrors(absl::Span< const double > coefficients, absl::Span< const double > lower_bounds, absl::Span< const double > upper_bounds, int64_t max_absolute_activity, double wanted_absolute_activity_precision, double *relative_coeff_error, double *scaled_sum_error)
bool RefIsPositive(int ref)
IntegerValue CeilRatio(IntegerValue dividend, IntegerValue positive_divisor)
bool ConvertMPModelProtoToCpModelProto(const SatParameters ¶ms, const MPModelProto &mp_model, CpModelProto *cp_model, SolverLogger *logger)
int64_t FindRationalFactor(double x, int64_t limit, double tolerance)
void ConvertBooleanProblemToLinearProgram(const LinearBooleanProblem &problem, glop::LinearProgram *lp)
Converts a Boolean optimization problem to its lp formulation.
bool MakeBoundsOfIntegerVariablesInteger(const SatParameters ¶ms, MPModelProto *mp_model, SolverLogger *logger)
void ChangeLargeBoundsToInfinity(double max_magnitude, MPModelProto *mp_model, SolverLogger *logger)
bool ScaleAndSetObjective(const SatParameters ¶ms, absl::Span< const std::pair< int, double > > objective, double objective_offset, bool maximize, CpModelProto *cp_model, SolverLogger *logger)
void ChangeOptimizationDirection(LinearBooleanProblem *problem)
std::vector< double > DetectImpliedIntegers(MPModelProto *mp_model, SolverLogger *logger)
void RemoveNearZeroTerms(const SatParameters ¶ms, MPModelProto *mp_model, SolverLogger *logger)
bool ConvertBinaryMPModelProtoToBooleanProblem(const MPModelProto &mp_model, LinearBooleanProblem *problem)
double ComputeTrueObjectiveLowerBound(const CpModelProto &model_proto_with_floating_point_objective, const CpObjectiveProto &integer_objective, const int64_t inner_integer_objective_lower_bound)
bool ConvertCpModelProtoToMPModelProto(const CpModelProto &input, MPModelProto *output)
constexpr Fractional kInfinity
Infinity for type Fractional.
std::vector< double > ScaleContinuousVariables(double scaling, double max_bound, MPModelProto *mp_model)
int NegatedRef(int ref)
Small utility functions to deal with negative variable/literal references.
bool MPModelProtoValidationBeforeConversion(const SatParameters ¶ms, const MPModelProto &mp_model, SolverLogger *logger)
In SWIG mode, we don't want anything besides these top-level includes.
int64_t CapAdd(int64_t x, int64_t y)
double GetBestScalingOfDoublesToInt64(absl::Span< const double > input, absl::Span< const double > lb, absl::Span< const double > ub, int64_t max_absolute_sum)
std::string ProtobufShortDebugString(const P &message)
int64_t CapProd(int64_t x, int64_t y)
void ComputeScalingErrors(absl::Span< const double > input, absl::Span< const double > lb, absl::Span< const double > ub, double scaling_factor, double *max_relative_coeff_error, double *max_scaled_sum_error)
std::string ProtobufDebugString(const P &message)
int64_t ComputeGcdOfRoundedDoubles(absl::Span< const double > x, double scaling_factor)
static int input(yyscan_t yyscanner)
#define SOLVER_LOG(logger,...)