1 // Copyright 2010-2021 Google LLC
2 // Licensed under the Apache License, Version 2.0 (the "License");
3 // you may not use this file except in compliance with the License.
4 // You may obtain a copy of the License at
5 //
6 // http://www.apache.org/licenses/LICENSE-2.0
7 //
8 // Unless required by applicable law or agreed to in writing, software
9 // distributed under the License is distributed on an "AS IS" BASIS,
10 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11 // See the License for the specific language governing permissions and
12 // limitations under the License.
13
14 #include "ortools/sat/lp_utils.h"
15
16 #include <stdlib.h>
17
18 #include <algorithm>
19 #include <cmath>
20 #include <cstdint>
21 #include <limits>
22 #include <string>
23 #include <vector>
24
25 #include "absl/strings/str_cat.h"
26 #include "ortools/base/int_type.h"
27 #include "ortools/base/integral_types.h"
28 #include "ortools/base/logging.h"
29 #include "ortools/glop/lp_solver.h"
30 #include "ortools/glop/parameters.pb.h"
31 #include "ortools/lp_data/lp_types.h"
32 #include "ortools/sat/boolean_problem.h"
33 #include "ortools/sat/cp_model_utils.h"
34 #include "ortools/sat/integer.h"
35 #include "ortools/sat/sat_base.h"
36 #include "ortools/util/fp_utils.h"
37
38 namespace operations_research {
39 namespace sat {
40
41 using glop::ColIndex;
42 using glop::Fractional;
43 using glop::kInfinity;
44 using glop::RowIndex;
45
46 using operations_research::MPConstraintProto;
47 using operations_research::MPModelProto;
48 using operations_research::MPVariableProto;
49
50 namespace {
51
ScaleConstraint(const std::vector<double> & var_scaling,MPConstraintProto * mp_constraint)52 void ScaleConstraint(const std::vector<double>& var_scaling,
53 MPConstraintProto* mp_constraint) {
54 const int num_terms = mp_constraint->coefficient_size();
55 for (int i = 0; i < num_terms; ++i) {
56 const int var_index = mp_constraint->var_index(i);
57 mp_constraint->set_coefficient(
58 i, mp_constraint->coefficient(i) / var_scaling[var_index]);
59 }
60 }
61
ApplyVarScaling(const std::vector<double> & var_scaling,MPModelProto * mp_model)62 void ApplyVarScaling(const std::vector<double>& var_scaling,
63 MPModelProto* mp_model) {
64 const int num_variables = mp_model->variable_size();
65 for (int i = 0; i < num_variables; ++i) {
66 const double scaling = var_scaling[i];
67 const MPVariableProto& mp_var = mp_model->variable(i);
68 const double old_lb = mp_var.lower_bound();
69 const double old_ub = mp_var.upper_bound();
70 const double old_obj = mp_var.objective_coefficient();
71 mp_model->mutable_variable(i)->set_lower_bound(old_lb * scaling);
72 mp_model->mutable_variable(i)->set_upper_bound(old_ub * scaling);
73 mp_model->mutable_variable(i)->set_objective_coefficient(old_obj / scaling);
74
75 // TODO(user): Make bounds of integer variable integer.
76 }
77 for (MPConstraintProto& mp_constraint : *mp_model->mutable_constraint()) {
78 ScaleConstraint(var_scaling, &mp_constraint);
79 }
80 for (MPGeneralConstraintProto& general_constraint :
81 *mp_model->mutable_general_constraint()) {
82 switch (general_constraint.general_constraint_case()) {
83 case MPGeneralConstraintProto::kIndicatorConstraint:
84 ScaleConstraint(var_scaling,
85 general_constraint.mutable_indicator_constraint()
86 ->mutable_constraint());
87 break;
88 case MPGeneralConstraintProto::kAndConstraint:
89 case MPGeneralConstraintProto::kOrConstraint:
90 // These constraints have only Boolean variables and no constants. They
91 // don't need scaling.
92 break;
93 default:
94 LOG(FATAL) << "Scaling unsupported for general constraint of type "
95 << general_constraint.general_constraint_case();
96 }
97 }
98 }
99
100 } // namespace
101
ScaleContinuousVariables(double scaling,double max_bound,MPModelProto * mp_model)102 std::vector<double> ScaleContinuousVariables(double scaling, double max_bound,
103 MPModelProto* mp_model) {
104 const int num_variables = mp_model->variable_size();
105 std::vector<double> var_scaling(num_variables, 1.0);
106 for (int i = 0; i < num_variables; ++i) {
107 if (mp_model->variable(i).is_integer()) continue;
108 const double lb = mp_model->variable(i).lower_bound();
109 const double ub = mp_model->variable(i).upper_bound();
110 const double magnitude = std::max(std::abs(lb), std::abs(ub));
111 if (magnitude == 0 || magnitude > max_bound) continue;
112 var_scaling[i] = std::min(scaling, max_bound / magnitude);
113 }
114 ApplyVarScaling(var_scaling, mp_model);
115 return var_scaling;
116 }
117
118 // This uses the best rational approximation of x via continuous fractions. It
119 // is probably not the best implementation, but according to the unit test, it
120 // seems to do the job.
FindRationalFactor(double x,int limit,double tolerance)121 int FindRationalFactor(double x, int limit, double tolerance) {
122 const double initial_x = x;
123 x = std::abs(x);
124 x -= std::floor(x);
125 int q = 1;
126 int prev_q = 0;
127 while (q < limit) {
128 if (std::abs(q * initial_x - std::round(q * initial_x)) < q * tolerance) {
129 return q;
130 }
131 x = 1 / x;
132 const int new_q = prev_q + static_cast<int>(std::floor(x)) * q;
133 prev_q = q;
134 q = new_q;
135 x -= std::floor(x);
136 }
137 return 0;
138 }
139
140 namespace {
141
142 // Returns a factor such that factor * var only need to take integer values to
143 // satisfy the given constraint. Return 0.0 if we didn't find such factor.
144 //
145 // Precondition: var must be the only non-integer in the given constraint.
GetIntegralityMultiplier(const MPModelProto & mp_model,const std::vector<double> & var_scaling,int var,int ct_index,double tolerance)146 double GetIntegralityMultiplier(const MPModelProto& mp_model,
147 const std::vector<double>& var_scaling, int var,
148 int ct_index, double tolerance) {
149 DCHECK(!mp_model.variable(var).is_integer());
150 const MPConstraintProto& ct = mp_model.constraint(ct_index);
151 double multiplier = 1.0;
152 double var_coeff = 0.0;
153 const double max_multiplier = 1e4;
154 for (int i = 0; i < ct.var_index().size(); ++i) {
155 if (var == ct.var_index(i)) {
156 var_coeff = ct.coefficient(i);
157 continue;
158 }
159
160 DCHECK(mp_model.variable(ct.var_index(i)).is_integer());
161 // This actually compute the smallest multiplier to make all other
162 // terms in the constraint integer.
163 const double coeff =
164 multiplier * ct.coefficient(i) / var_scaling[ct.var_index(i)];
165 multiplier *=
166 FindRationalFactor(coeff, /*limit=*/100, multiplier * tolerance);
167 if (multiplier == 0 || multiplier > max_multiplier) return 0.0;
168 }
169 DCHECK_NE(var_coeff, 0.0);
170
171 // The constraint bound need to be infinite or integer.
172 for (const double bound : {ct.lower_bound(), ct.upper_bound()}) {
173 if (!std::isfinite(bound)) continue;
174 if (std::abs(std::round(bound * multiplier) - bound * multiplier) >
175 tolerance * multiplier) {
176 return 0.0;
177 }
178 }
179 return std::abs(multiplier * var_coeff);
180 }
181
182 } // namespace
183
MakeBoundsOfIntegerVariablesInteger(const SatParameters & params,MPModelProto * mp_model,SolverLogger * logger)184 bool MakeBoundsOfIntegerVariablesInteger(const SatParameters& params,
185 MPModelProto* mp_model,
186 SolverLogger* logger) {
187 const int num_variables = mp_model->variable_size();
188 const double tolerance = params.mip_wanted_precision();
189 int64_t num_changes = 0;
190 for (int i = 0; i < num_variables; ++i) {
191 const MPVariableProto& mp_var = mp_model->variable(i);
192 if (!mp_var.is_integer()) continue;
193
194 const double lb = mp_var.lower_bound();
195 const double new_lb = std::isfinite(lb) ? std::ceil(lb - tolerance) : lb;
196 if (lb != new_lb) {
197 ++num_changes;
198 mp_model->mutable_variable(i)->set_lower_bound(new_lb);
199 }
200
201 const double ub = mp_var.upper_bound();
202 const double new_ub = std::isfinite(ub) ? std::floor(ub + tolerance) : ub;
203 if (ub != new_ub) {
204 ++num_changes;
205 mp_model->mutable_variable(i)->set_upper_bound(new_ub);
206 }
207
208 if (new_ub < new_lb) {
209 SOLVER_LOG(logger, "Empty domain for integer variable #", i, ": [", lb,
210 ",", ub, "]");
211 return false;
212 }
213 }
214
215 if (num_changes > 0) {
216 SOLVER_LOG(logger, "Changed ", num_changes,
217 " bounds of integer variables to integer values");
218 }
219 return true;
220 }
221
RemoveNearZeroTerms(const SatParameters & params,MPModelProto * mp_model,SolverLogger * logger)222 void RemoveNearZeroTerms(const SatParameters& params, MPModelProto* mp_model,
223 SolverLogger* logger) {
224 const int num_variables = mp_model->variable_size();
225
226 // Compute for each variable its current maximum magnitude. Note that we will
227 // only scale variable with a coefficient >= 1, so it is safe to use this
228 // bound.
229 std::vector<double> max_bounds(num_variables);
230 for (int i = 0; i < num_variables; ++i) {
231 double value = std::abs(mp_model->variable(i).lower_bound());
232 value = std::max(value, std::abs(mp_model->variable(i).upper_bound()));
233 value = std::min(value, params.mip_max_bound());
234 max_bounds[i] = value;
235 }
236
237 // We want the maximum absolute error while setting coefficients to zero to
238 // not exceed our mip wanted precision. So for a binary variable we might set
239 // to zero coefficient around 1e-7. But for large domain, we need lower coeff
240 // than that, around 1e-12 with the default params.mip_max_bound(). This also
241 // depends on the size of the constraint.
242 int64_t num_removed = 0;
243 double largest_removed = 0.0;
244 const int num_constraints = mp_model->constraint_size();
245 for (int c = 0; c < num_constraints; ++c) {
246 MPConstraintProto* ct = mp_model->mutable_constraint(c);
247 int new_size = 0;
248 const int size = ct->var_index().size();
249 if (size == 0) continue;
250 const double threshold =
251 params.mip_wanted_precision() / static_cast<double>(size);
252 for (int i = 0; i < size; ++i) {
253 const int var = ct->var_index(i);
254 const double coeff = ct->coefficient(i);
255 if (std::abs(coeff) * max_bounds[var] < threshold) {
256 largest_removed = std::max(largest_removed, std::abs(coeff));
257 continue;
258 }
259 ct->set_var_index(new_size, var);
260 ct->set_coefficient(new_size, coeff);
261 ++new_size;
262 }
263 num_removed += size - new_size;
264 ct->mutable_var_index()->Truncate(new_size);
265 ct->mutable_coefficient()->Truncate(new_size);
266 }
267
268 if (num_removed > 0) {
269 SOLVER_LOG(logger, "Removed ", num_removed,
270 " near zero terms with largest magnitude of ", largest_removed,
271 ".");
272 }
273 }
274
MPModelProtoValidationBeforeConversion(const SatParameters & params,const MPModelProto & mp_model,SolverLogger * logger)275 bool MPModelProtoValidationBeforeConversion(const SatParameters& params,
276 const MPModelProto& mp_model,
277 SolverLogger* logger) {
278 // Abort if there is constraint type we don't currently support.
279 for (const MPGeneralConstraintProto& general_constraint :
280 mp_model.general_constraint()) {
281 switch (general_constraint.general_constraint_case()) {
282 case MPGeneralConstraintProto::kIndicatorConstraint:
283 break;
284 case MPGeneralConstraintProto::kAndConstraint:
285 break;
286 case MPGeneralConstraintProto::kOrConstraint:
287 break;
288 default:
289 SOLVER_LOG(logger, "General constraints of type ",
290 general_constraint.general_constraint_case(),
291 " are not supported.");
292 return false;
293 }
294 }
295
296 // Abort if finite variable bounds or objective is too large.
297 const double threshold = params.mip_max_valid_magnitude();
298 const int num_variables = mp_model.variable_size();
299 for (int i = 0; i < num_variables; ++i) {
300 const MPVariableProto& var = mp_model.variable(i);
301 if ((std::isfinite(var.lower_bound()) &&
302 std::abs(var.lower_bound()) > threshold) ||
303 (std::isfinite(var.upper_bound()) &&
304 std::abs(var.upper_bound()) > threshold)) {
305 SOLVER_LOG(logger, "Variable bounds are too large [", var.lower_bound(),
306 ",", var.upper_bound(), "]");
307 return false;
308 }
309 if (std::abs(var.objective_coefficient()) > threshold) {
310 SOLVER_LOG(logger, "Objective coefficient is too large: ",
311 var.objective_coefficient());
312 return false;
313 }
314 }
315
316 // Abort if finite constraint bounds or coefficients are too large.
317 const int num_constraints = mp_model.constraint_size();
318 for (int c = 0; c < num_constraints; ++c) {
319 const MPConstraintProto& ct = mp_model.constraint(c);
320 if ((std::isfinite(ct.lower_bound()) &&
321 std::abs(ct.lower_bound()) > threshold) ||
322 (std::isfinite(ct.upper_bound()) &&
323 std::abs(ct.upper_bound()) > threshold)) {
324 SOLVER_LOG(logger, "Constraint bounds are too large [", ct.lower_bound(),
325 ",", ct.upper_bound(), "]");
326 return false;
327 }
328 for (const double coeff : ct.coefficient()) {
329 if (std::abs(coeff) > threshold) {
330 SOLVER_LOG(logger, "Constraint coefficient is too large: ", coeff);
331 return false;
332 }
333 }
334 }
335
336 return true;
337 }
338
DetectImpliedIntegers(MPModelProto * mp_model,SolverLogger * logger)339 std::vector<double> DetectImpliedIntegers(MPModelProto* mp_model,
340 SolverLogger* logger) {
341 const int num_variables = mp_model->variable_size();
342 std::vector<double> var_scaling(num_variables, 1.0);
343
344 int initial_num_integers = 0;
345 for (int i = 0; i < num_variables; ++i) {
346 if (mp_model->variable(i).is_integer()) ++initial_num_integers;
347 }
348 VLOG(1) << "Initial num integers: " << initial_num_integers;
349
350 // We will process all equality constraints with exactly one non-integer.
351 const double tolerance = 1e-6;
352 std::vector<int> constraint_queue;
353
354 const int num_constraints = mp_model->constraint_size();
355 std::vector<int> constraint_to_num_non_integer(num_constraints, 0);
356 std::vector<std::vector<int>> var_to_constraints(num_variables);
357 for (int i = 0; i < num_constraints; ++i) {
358 const MPConstraintProto& mp_constraint = mp_model->constraint(i);
359
360 for (const int var : mp_constraint.var_index()) {
361 if (!mp_model->variable(var).is_integer()) {
362 var_to_constraints[var].push_back(i);
363 constraint_to_num_non_integer[i]++;
364 }
365 }
366 if (constraint_to_num_non_integer[i] == 1) {
367 constraint_queue.push_back(i);
368 }
369 }
370 VLOG(1) << "Initial constraint queue: " << constraint_queue.size() << " / "
371 << num_constraints;
372
373 int num_detected = 0;
374 double max_scaling = 0.0;
375 auto scale_and_mark_as_integer = [&](int var, double scaling) mutable {
376 CHECK_NE(var, -1);
377 CHECK(!mp_model->variable(var).is_integer());
378 CHECK_EQ(var_scaling[var], 1.0);
379 if (scaling != 1.0) {
380 VLOG(2) << "Scaled " << var << " by " << scaling;
381 }
382
383 ++num_detected;
384 max_scaling = std::max(max_scaling, scaling);
385
386 // Scale the variable right away and mark it as implied integer.
387 // Note that the constraints will be scaled later.
388 var_scaling[var] = scaling;
389 mp_model->mutable_variable(var)->set_is_integer(true);
390
391 // Update the queue of constraints with a single non-integer.
392 for (const int ct_index : var_to_constraints[var]) {
393 constraint_to_num_non_integer[ct_index]--;
394 if (constraint_to_num_non_integer[ct_index] == 1) {
395 constraint_queue.push_back(ct_index);
396 }
397 }
398 };
399
400 int num_fail_due_to_rhs = 0;
401 int num_fail_due_to_large_multiplier = 0;
402 int num_processed_constraints = 0;
403 while (!constraint_queue.empty()) {
404 const int top_ct_index = constraint_queue.back();
405 constraint_queue.pop_back();
406
407 // The non integer variable was already made integer by one other
408 // constraint.
409 if (constraint_to_num_non_integer[top_ct_index] == 0) continue;
410
411 // Ignore non-equality here.
412 const MPConstraintProto& ct = mp_model->constraint(top_ct_index);
413 if (ct.lower_bound() + tolerance < ct.upper_bound()) continue;
414
415 ++num_processed_constraints;
416
417 // This will be set to the unique non-integer term of this constraint.
418 int var = -1;
419 double var_coeff;
420
421 // We are looking for a "multiplier" so that the unique non-integer term
422 // in this constraint (i.e. var * var_coeff) times this multiplier is an
423 // integer.
424 //
425 // If this is set to zero or becomes too large, we fail to detect a new
426 // implied integer and ignore this constraint.
427 double multiplier = 1.0;
428 const double max_multiplier = 1e4;
429
430 for (int i = 0; i < ct.var_index().size(); ++i) {
431 if (!mp_model->variable(ct.var_index(i)).is_integer()) {
432 CHECK_EQ(var, -1);
433 var = ct.var_index(i);
434 var_coeff = ct.coefficient(i);
435 } else {
436 // This actually compute the smallest multiplier to make all other
437 // terms in the constraint integer.
438 const double coeff =
439 multiplier * ct.coefficient(i) / var_scaling[ct.var_index(i)];
440 multiplier *=
441 FindRationalFactor(coeff, /*limit=*/100, multiplier * tolerance);
442 if (multiplier == 0 || multiplier > max_multiplier) {
443 break;
444 }
445 }
446 }
447
448 if (multiplier == 0 || multiplier > max_multiplier) {
449 ++num_fail_due_to_large_multiplier;
450 continue;
451 }
452
453 // These "rhs" fail could be handled by shifting the variable.
454 const double rhs = ct.lower_bound();
455 if (std::abs(std::round(rhs * multiplier) - rhs * multiplier) >
456 tolerance * multiplier) {
457 ++num_fail_due_to_rhs;
458 continue;
459 }
460
461 // We want to multiply the variable so that it is integer. We know that
462 // coeff * multiplier is an integer, so we just multiply by that.
463 //
464 // But if a variable appear in more than one equality, we want to find the
465 // smallest integrality factor! See diameterc-msts-v40a100d5i.mps
466 // for an instance of this.
467 double best_scaling = std::abs(var_coeff * multiplier);
468 for (const int ct_index : var_to_constraints[var]) {
469 if (ct_index == top_ct_index) continue;
470 if (constraint_to_num_non_integer[ct_index] != 1) continue;
471
472 // Ignore non-equality here.
473 const MPConstraintProto& ct = mp_model->constraint(top_ct_index);
474 if (ct.lower_bound() + tolerance < ct.upper_bound()) continue;
475
476 const double multiplier = GetIntegralityMultiplier(
477 *mp_model, var_scaling, var, ct_index, tolerance);
478 if (multiplier != 0.0 && multiplier < best_scaling) {
479 best_scaling = multiplier;
480 }
481 }
482
483 scale_and_mark_as_integer(var, best_scaling);
484 }
485
486 // Process continuous variables that only appear as the unique non integer
487 // in a set of non-equality constraints.
488 //
489 // Note that turning to integer such variable cannot in turn trigger new
490 // integer detection, so there is no point doing that in a loop.
491 int num_in_inequalities = 0;
492 int num_to_be_handled = 0;
493 for (int var = 0; var < num_variables; ++var) {
494 if (mp_model->variable(var).is_integer()) continue;
495
496 // This should be presolved and not happen.
497 if (var_to_constraints[var].empty()) continue;
498
499 bool ok = true;
500 for (const int ct_index : var_to_constraints[var]) {
501 if (constraint_to_num_non_integer[ct_index] != 1) {
502 ok = false;
503 break;
504 }
505 }
506 if (!ok) continue;
507
508 std::vector<double> scaled_coeffs;
509 for (const int ct_index : var_to_constraints[var]) {
510 const double multiplier = GetIntegralityMultiplier(
511 *mp_model, var_scaling, var, ct_index, tolerance);
512 if (multiplier == 0.0) {
513 ok = false;
514 break;
515 }
516 scaled_coeffs.push_back(multiplier);
517 }
518 if (!ok) continue;
519
520 // The situation is a bit tricky here, we have a bunch of coeffs c_i, and we
521 // know that X * c_i can take integer value without changing the constraint
522 // i meaning.
523 //
524 // For now we take the min, and scale only if all c_i / min are integer.
525 double scaling = scaled_coeffs[0];
526 for (const double c : scaled_coeffs) {
527 scaling = std::min(scaling, c);
528 }
529 CHECK_GT(scaling, 0.0);
530 for (const double c : scaled_coeffs) {
531 const double fraction = c / scaling;
532 if (std::abs(std::round(fraction) - fraction) > tolerance) {
533 ok = false;
534 break;
535 }
536 }
537 if (!ok) {
538 // TODO(user): be smarter! we should be able to handle these cases.
539 ++num_to_be_handled;
540 continue;
541 }
542
543 // Tricky, we also need the bound of the scaled variable to be integer.
544 for (const double bound : {mp_model->variable(var).lower_bound(),
545 mp_model->variable(var).upper_bound()}) {
546 if (!std::isfinite(bound)) continue;
547 if (std::abs(std::round(bound * scaling) - bound * scaling) >
548 tolerance * scaling) {
549 ok = false;
550 break;
551 }
552 }
553 if (!ok) {
554 // TODO(user): If we scale more we migth be able to turn it into an
555 // integer.
556 ++num_to_be_handled;
557 continue;
558 }
559
560 ++num_in_inequalities;
561 scale_and_mark_as_integer(var, scaling);
562 }
563 VLOG(1) << "num_new_integer: " << num_detected
564 << " num_processed_constraints: " << num_processed_constraints
565 << " num_rhs_fail: " << num_fail_due_to_rhs
566 << " num_multiplier_fail: " << num_fail_due_to_large_multiplier;
567
568 if (num_to_be_handled > 0) {
569 SOLVER_LOG(logger, "Missed ", num_to_be_handled,
570 " potential implied integer.");
571 }
572
573 const int num_integers = initial_num_integers + num_detected;
574 SOLVER_LOG(logger, "Num integers: ", num_integers, "/", num_variables,
575 " (implied: ", num_detected,
576 " in_inequalities: ", num_in_inequalities,
577 " max_scaling: ", max_scaling, ")",
578 (num_integers == num_variables ? " [IP] " : " [MIP] "));
579
580 ApplyVarScaling(var_scaling, mp_model);
581 return var_scaling;
582 }
583
584 namespace {
585
586 // We use a class to reuse the temporary memory.
587 struct ConstraintScaler {
588 // Scales an individual constraint.
589 ConstraintProto* AddConstraint(const MPModelProto& mp_model,
590 const MPConstraintProto& mp_constraint,
591 CpModelProto* cp_model);
592
593 double max_relative_coeff_error = 0.0;
594 double max_relative_rhs_error = 0.0;
595 double max_scaling_factor = 0.0;
596
597 double wanted_precision = 1e-6;
598 int64_t scaling_target = int64_t{1} << 50;
599 std::vector<int> var_indices;
600 std::vector<double> coefficients;
601 std::vector<double> lower_bounds;
602 std::vector<double> upper_bounds;
603 };
604
605 namespace {
606
FindFractionalScaling(const std::vector<double> & coefficients,double tolerance)607 double FindFractionalScaling(const std::vector<double>& coefficients,
608 double tolerance) {
609 double multiplier = 1.0;
610 for (const double coeff : coefficients) {
611 multiplier *=
612 FindRationalFactor(coeff, /*limit=*/1e8, multiplier * tolerance);
613 if (multiplier == 0.0) break;
614 }
615 return multiplier;
616 }
617
618 } // namespace
619
AddConstraint(const MPModelProto & mp_model,const MPConstraintProto & mp_constraint,CpModelProto * cp_model)620 ConstraintProto* ConstraintScaler::AddConstraint(
621 const MPModelProto& mp_model, const MPConstraintProto& mp_constraint,
622 CpModelProto* cp_model) {
623 if (mp_constraint.lower_bound() == -kInfinity &&
624 mp_constraint.upper_bound() == kInfinity) {
625 return nullptr;
626 }
627
628 auto* constraint = cp_model->add_constraints();
629 constraint->set_name(mp_constraint.name());
630 auto* arg = constraint->mutable_linear();
631
632 // First scale the coefficients of the constraints so that the constraint
633 // sum can always be computed without integer overflow.
634 var_indices.clear();
635 coefficients.clear();
636 lower_bounds.clear();
637 upper_bounds.clear();
638 const int num_coeffs = mp_constraint.coefficient_size();
639 for (int i = 0; i < num_coeffs; ++i) {
640 const auto& var_proto = cp_model->variables(mp_constraint.var_index(i));
641 const int64_t lb = var_proto.domain(0);
642 const int64_t ub = var_proto.domain(var_proto.domain_size() - 1);
643 if (lb == 0 && ub == 0) continue;
644
645 const double coeff = mp_constraint.coefficient(i);
646 if (coeff == 0.0) continue;
647
648 var_indices.push_back(mp_constraint.var_index(i));
649 coefficients.push_back(coeff);
650 lower_bounds.push_back(lb);
651 upper_bounds.push_back(ub);
652 }
653
654 // We compute the worst case error relative to the magnitude of the bounds.
655 Fractional lb = mp_constraint.lower_bound();
656 Fractional ub = mp_constraint.upper_bound();
657 const double ct_norm = std::max(1.0, std::min(std::abs(lb), std::abs(ub)));
658
659 double scaling_factor = GetBestScalingOfDoublesToInt64(
660 coefficients, lower_bounds, upper_bounds, scaling_target);
661 if (scaling_factor == 0.0) {
662 // TODO(user): Report error properly instead of ignoring constraint. Note
663 // however that this likely indicate a coefficient of inf in the constraint,
664 // so we should probably abort before reaching here.
665 LOG(DFATAL) << "Scaling factor of zero while scaling constraint: "
666 << mp_constraint.ShortDebugString();
667 return nullptr;
668 }
669
670 // Returns the smallest factor of the form 2^i that gives us a relative sum
671 // error of wanted_precision and still make sure we will have no integer
672 // overflow.
673 //
674 // TODO(user): Make this faster.
675 double x = std::min(scaling_factor, 1.0);
676 double relative_coeff_error;
677 double scaled_sum_error;
678 for (; x <= scaling_factor; x *= 2) {
679 ComputeScalingErrors(coefficients, lower_bounds, upper_bounds, x,
680 &relative_coeff_error, &scaled_sum_error);
681 if (scaled_sum_error < wanted_precision * x * ct_norm) break;
682 }
683 scaling_factor = x;
684
685 // Because we deal with an approximate input, scaling with a power of 2 might
686 // not be the best choice. It is also possible user used rational coeff and
687 // then converted them to double (1/2, 1/3, 4/5, etc...). This scaling will
688 // recover such rational input and might result in a smaller overall
689 // coefficient which is good.
690 const double integer_factor = FindFractionalScaling(coefficients, 1e-8);
691 if (integer_factor != 0 && integer_factor < scaling_factor) {
692 ComputeScalingErrors(coefficients, lower_bounds, upper_bounds, x,
693 &relative_coeff_error, &scaled_sum_error);
694 if (scaled_sum_error < wanted_precision * integer_factor * ct_norm) {
695 scaling_factor = integer_factor;
696 }
697 }
698
699 const int64_t gcd = ComputeGcdOfRoundedDoubles(coefficients, scaling_factor);
700 max_relative_coeff_error =
701 std::max(relative_coeff_error, max_relative_coeff_error);
702 max_scaling_factor = std::max(scaling_factor / gcd, max_scaling_factor);
703
704 // We do not relax the constraint bound if all variables are integer and
705 // we made no error at all during our scaling.
706 bool relax_bound = scaled_sum_error > 0;
707
708 for (int i = 0; i < coefficients.size(); ++i) {
709 const double scaled_value = coefficients[i] * scaling_factor;
710 const int64_t value = static_cast<int64_t>(std::round(scaled_value)) / gcd;
711 if (value != 0) {
712 if (!mp_model.variable(var_indices[i]).is_integer()) {
713 relax_bound = true;
714 }
715 arg->add_vars(var_indices[i]);
716 arg->add_coeffs(value);
717 }
718 }
719 max_relative_rhs_error = std::max(
720 max_relative_rhs_error, scaled_sum_error / (scaling_factor * ct_norm));
721
722 // Add the constraint bounds. Because we are sure the scaled constraint fit
723 // on an int64_t, if the scaled bounds are too large, the constraint is either
724 // always true or always false.
725 if (relax_bound) {
726 lb -= std::max(std::abs(lb), 1.0) * wanted_precision;
727 }
728 const Fractional scaled_lb = std::ceil(lb * scaling_factor);
729 if (lb == kInfinity || scaled_lb >= std::numeric_limits<int64_t>::max()) {
730 // Corner case: infeasible model.
731 arg->add_domain(std::numeric_limits<int64_t>::max());
732 } else if (lb == -kInfinity ||
733 scaled_lb <= std::numeric_limits<int64_t>::min()) {
734 arg->add_domain(std::numeric_limits<int64_t>::min());
735 } else {
736 arg->add_domain(CeilRatio(IntegerValue(static_cast<int64_t>(scaled_lb)),
737 IntegerValue(gcd))
738 .value());
739 }
740
741 if (relax_bound) {
742 ub += std::max(std::abs(ub), 1.0) * wanted_precision;
743 }
744 const Fractional scaled_ub = std::floor(ub * scaling_factor);
745 if (ub == -kInfinity || scaled_ub <= std::numeric_limits<int64_t>::min()) {
746 // Corner case: infeasible model.
747 arg->add_domain(std::numeric_limits<int64_t>::min());
748 } else if (ub == kInfinity ||
749 scaled_ub >= std::numeric_limits<int64_t>::max()) {
750 arg->add_domain(std::numeric_limits<int64_t>::max());
751 } else {
752 arg->add_domain(FloorRatio(IntegerValue(static_cast<int64_t>(scaled_ub)),
753 IntegerValue(gcd))
754 .value());
755 }
756
757 return constraint;
758 }
759
760 } // namespace
761
ConvertMPModelProtoToCpModelProto(const SatParameters & params,const MPModelProto & mp_model,CpModelProto * cp_model,SolverLogger * logger)762 bool ConvertMPModelProtoToCpModelProto(const SatParameters& params,
763 const MPModelProto& mp_model,
764 CpModelProto* cp_model,
765 SolverLogger* logger) {
766 CHECK(cp_model != nullptr);
767 cp_model->Clear();
768 cp_model->set_name(mp_model.name());
769
770 // To make sure we cannot have integer overflow, we use this bound for any
771 // unbounded variable.
772 //
773 // TODO(user): This could be made larger if needed, so be smarter if we have
774 // MIP problem that we cannot "convert" because of this. Note however than we
775 // cannot go that much further because we need to make sure we will not run
776 // into overflow if we add a big linear combination of such variables. It
777 // should always be possible for a user to scale its problem so that all
778 // relevant quantities are a couple of millions. A LP/MIP solver have a
779 // similar condition in disguise because problem with a difference of more
780 // than 6 magnitudes between the variable values will likely run into numeric
781 // trouble.
782 const int64_t kMaxVariableBound =
783 static_cast<int64_t>(params.mip_max_bound());
784
785 int num_truncated_bounds = 0;
786 int num_small_domains = 0;
787 const int64_t kSmallDomainSize = 1000;
788 const double kWantedPrecision = params.mip_wanted_precision();
789
790 // Add the variables.
791 const int num_variables = mp_model.variable_size();
792 for (int i = 0; i < num_variables; ++i) {
793 const MPVariableProto& mp_var = mp_model.variable(i);
794 IntegerVariableProto* cp_var = cp_model->add_variables();
795 cp_var->set_name(mp_var.name());
796
797 // Deal with the corner case of a domain far away from zero.
798 //
799 // TODO(user): We could avoid these cases by shifting the domain of
800 // all variables to contain zero. This should also lead to a better scaling,
801 // but it has some complications with integer variables and require some
802 // post-solve.
803 if (mp_var.lower_bound() > static_cast<double>(kMaxVariableBound) ||
804 mp_var.upper_bound() < static_cast<double>(-kMaxVariableBound)) {
805 SOLVER_LOG(logger, "Error: variable ", mp_var.DebugString(),
806 " is outside [-mip_max_bound..mip_max_bound]");
807 return false;
808 }
809
810 // Note that we must process the lower bound first.
811 for (const bool lower : {true, false}) {
812 const double bound = lower ? mp_var.lower_bound() : mp_var.upper_bound();
813 if (std::abs(bound) >= static_cast<double>(kMaxVariableBound)) {
814 ++num_truncated_bounds;
815 cp_var->add_domain(bound < 0 ? -kMaxVariableBound : kMaxVariableBound);
816 continue;
817 }
818
819 // Note that the cast is "perfect" because we forbid large values.
820 cp_var->add_domain(
821 static_cast<int64_t>(lower ? std::ceil(bound - kWantedPrecision)
822 : std::floor(bound + kWantedPrecision)));
823 }
824
825 if (cp_var->domain(0) > cp_var->domain(1)) {
826 LOG(WARNING) << "Variable #" << i << " cannot take integer value. "
827 << mp_var.ShortDebugString();
828 return false;
829 }
830
831 // Notify if a continuous variable has a small domain as this is likely to
832 // make an all integer solution far from a continuous one.
833 if (!mp_var.is_integer() && cp_var->domain(0) != cp_var->domain(1) &&
834 cp_var->domain(1) - cp_var->domain(0) < kSmallDomainSize) {
835 ++num_small_domains;
836 }
837 }
838
839 if (num_truncated_bounds > 0) {
840 SOLVER_LOG(logger, "Warning: ", num_truncated_bounds,
841 " bounds were truncated to ", kMaxVariableBound, ".");
842 }
843 if (num_small_domains > 0) {
844 SOLVER_LOG(logger, "Warning: ", num_small_domains,
845 " continuous variable domain with fewer than ", kSmallDomainSize,
846 " values.");
847 }
848
849 ConstraintScaler scaler;
850 const int64_t kScalingTarget = int64_t{1}
851 << params.mip_max_activity_exponent();
852 scaler.wanted_precision = kWantedPrecision;
853 scaler.scaling_target = kScalingTarget;
854
855 // Add the constraints. We scale each of them individually.
856 for (const MPConstraintProto& mp_constraint : mp_model.constraint()) {
857 scaler.AddConstraint(mp_model, mp_constraint, cp_model);
858 }
859 for (const MPGeneralConstraintProto& general_constraint :
860 mp_model.general_constraint()) {
861 switch (general_constraint.general_constraint_case()) {
862 case MPGeneralConstraintProto::kIndicatorConstraint: {
863 const auto& indicator_constraint =
864 general_constraint.indicator_constraint();
865 const MPConstraintProto& mp_constraint =
866 indicator_constraint.constraint();
867 ConstraintProto* ct =
868 scaler.AddConstraint(mp_model, mp_constraint, cp_model);
869 if (ct == nullptr) continue;
870
871 // Add the indicator.
872 const int var = indicator_constraint.var_index();
873 const int value = indicator_constraint.var_value();
874 ct->add_enforcement_literal(value == 1 ? var : NegatedRef(var));
875 break;
876 }
877 case MPGeneralConstraintProto::kAndConstraint: {
878 const auto& and_constraint = general_constraint.and_constraint();
879 const std::string& name = general_constraint.name();
880
881 ConstraintProto* ct_pos = cp_model->add_constraints();
882 ct_pos->set_name(name.empty() ? "" : absl::StrCat(name, "_pos"));
883 ct_pos->add_enforcement_literal(and_constraint.resultant_var_index());
884 *ct_pos->mutable_bool_and()->mutable_literals() =
885 and_constraint.var_index();
886
887 ConstraintProto* ct_neg = cp_model->add_constraints();
888 ct_neg->set_name(name.empty() ? "" : absl::StrCat(name, "_neg"));
889 ct_neg->add_enforcement_literal(
890 NegatedRef(and_constraint.resultant_var_index()));
891 for (const int var_index : and_constraint.var_index()) {
892 ct_neg->mutable_bool_or()->add_literals(NegatedRef(var_index));
893 }
894 break;
895 }
896 case MPGeneralConstraintProto::kOrConstraint: {
897 const auto& or_constraint = general_constraint.or_constraint();
898 const std::string& name = general_constraint.name();
899
900 ConstraintProto* ct_pos = cp_model->add_constraints();
901 ct_pos->set_name(name.empty() ? "" : absl::StrCat(name, "_pos"));
902 ct_pos->add_enforcement_literal(or_constraint.resultant_var_index());
903 *ct_pos->mutable_bool_or()->mutable_literals() =
904 or_constraint.var_index();
905
906 ConstraintProto* ct_neg = cp_model->add_constraints();
907 ct_neg->set_name(name.empty() ? "" : absl::StrCat(name, "_neg"));
908 ct_neg->add_enforcement_literal(
909 NegatedRef(or_constraint.resultant_var_index()));
910 for (const int var_index : or_constraint.var_index()) {
911 ct_neg->mutable_bool_and()->add_literals(NegatedRef(var_index));
912 }
913 break;
914 }
915 default:
916 LOG(ERROR) << "Can't convert general constraints of type "
917 << general_constraint.general_constraint_case()
918 << " to CpModelProto.";
919 return false;
920 }
921 }
922
923 // Display the error/scaling on the constraints.
924 SOLVER_LOG(logger, "Maximum constraint coefficient relative error: ",
925 scaler.max_relative_coeff_error);
926 SOLVER_LOG(logger, "Maximum constraint worst-case activity relative error: ",
927 scaler.max_relative_rhs_error,
928 (scaler.max_relative_rhs_error > params.mip_check_precision()
929 ? " [Potentially IMPRECISE]"
930 : ""));
931 SOLVER_LOG(logger,
932 "Maximum constraint scaling factor: ", scaler.max_scaling_factor);
933
934 // Since cp_model support a floating point objective, we use that. This will
935 // allow us to scale the objective a bit later so we can potentially do more
936 // domain reduction first.
937 auto* float_objective = cp_model->mutable_floating_point_objective();
938 float_objective->set_maximize(mp_model.maximize());
939 float_objective->set_offset(mp_model.objective_offset());
940 for (int i = 0; i < num_variables; ++i) {
941 const MPVariableProto& mp_var = mp_model.variable(i);
942 if (mp_var.objective_coefficient() != 0.0) {
943 float_objective->add_vars(i);
944 float_objective->add_coeffs(mp_var.objective_coefficient());
945 }
946 }
947
948 // If the objective is fixed to zero, we consider there is none.
949 if (float_objective->offset() == 0 && float_objective->vars().empty()) {
950 cp_model->clear_floating_point_objective();
951 }
952 return true;
953 }
954
ScaleAndSetObjective(const SatParameters & params,const std::vector<std::pair<int,double>> & objective,double objective_offset,bool maximize,CpModelProto * cp_model,SolverLogger * logger)955 bool ScaleAndSetObjective(const SatParameters& params,
956 const std::vector<std::pair<int, double>>& objective,
957 double objective_offset, bool maximize,
958 CpModelProto* cp_model, SolverLogger* logger) {
959 // Make sure the objective is currently empty.
960 cp_model->clear_objective();
961
962 // We filter constant terms and compute some needed quantities.
963 std::vector<int> var_indices;
964 std::vector<double> coefficients;
965 std::vector<double> lower_bounds;
966 std::vector<double> upper_bounds;
967 double min_magnitude = std::numeric_limits<double>::infinity();
968 double max_magnitude = 0.0;
969 double l1_norm = 0.0;
970 for (const auto [var, coeff] : objective) {
971 const auto& var_proto = cp_model->variables(var);
972 const int64_t lb = var_proto.domain(0);
973 const int64_t ub = var_proto.domain(var_proto.domain_size() - 1);
974 if (lb == ub) {
975 if (lb != 0) objective_offset += lb * coeff;
976 continue;
977 }
978 var_indices.push_back(var);
979 coefficients.push_back(coeff);
980 lower_bounds.push_back(lb);
981 upper_bounds.push_back(ub);
982
983 min_magnitude = std::min(min_magnitude, std::abs(coeff));
984 max_magnitude = std::max(max_magnitude, std::abs(coeff));
985 l1_norm += std::abs(coeff);
986 }
987
988 if (coefficients.empty() && objective_offset == 0.0) return true;
989
990 if (!coefficients.empty()) {
991 const double average_magnitude =
992 l1_norm / static_cast<double>(coefficients.size());
993 SOLVER_LOG(logger, "[Scaling] Floating point objective has ",
994 coefficients.size(), " terms with magnitude in [", min_magnitude,
995 ", ", max_magnitude, "] average = ", average_magnitude);
996 }
997
998 // Start by computing the largest possible factor that gives something that
999 // will not cause overflow when evaluating the objective value.
1000 const int64_t kScalingTarget = int64_t{1}
1001 << params.mip_max_activity_exponent();
1002 double scaling_factor = GetBestScalingOfDoublesToInt64(
1003 coefficients, lower_bounds, upper_bounds, kScalingTarget);
1004 if (scaling_factor == 0.0) {
1005 LOG(ERROR) << "Scaling factor of zero while scaling objective! This "
1006 "likely indicate an infinite coefficient in the objective.";
1007 return false;
1008 }
1009
1010 // Returns the smallest factor of the form 2^i that gives us an absolute
1011 // error of kWantedPrecision and still make sure we will have no integer
1012 // overflow.
1013 //
1014 // TODO(user): Make this faster.
1015 double x = std::min(scaling_factor, 1.0);
1016 double relative_coeff_error;
1017 double scaled_sum_error;
1018 const double kWantedPrecision =
1019 std::max(params.mip_wanted_precision(), params.absolute_gap_limit());
1020 for (; x <= scaling_factor; x *= 2) {
1021 ComputeScalingErrors(coefficients, lower_bounds, upper_bounds, x,
1022 &relative_coeff_error, &scaled_sum_error);
1023 if (scaled_sum_error < kWantedPrecision * x) break;
1024 }
1025 scaling_factor = x;
1026
1027 // Same remark as for the constraint scaling code.
1028 // TODO(user): Extract common code.
1029 const double integer_factor = FindFractionalScaling(coefficients, 1e-8);
1030 if (integer_factor != 0 && integer_factor < scaling_factor) {
1031 ComputeScalingErrors(coefficients, lower_bounds, upper_bounds, x,
1032 &relative_coeff_error, &scaled_sum_error);
1033 if (scaled_sum_error < kWantedPrecision * integer_factor) {
1034 scaling_factor = integer_factor;
1035 }
1036 }
1037
1038 const int64_t gcd = ComputeGcdOfRoundedDoubles(coefficients, scaling_factor);
1039
1040 // Display the objective error/scaling.
1041 SOLVER_LOG(logger, "[Scaling] Objective coefficient relative error: ",
1042 relative_coeff_error);
1043 SOLVER_LOG(logger, "[Scaling] Objective worst-case absolute error: ",
1044 scaled_sum_error / scaling_factor);
1045 SOLVER_LOG(logger,
1046 "[Scaling] Objective scaling factor: ", scaling_factor / gcd);
1047
1048 // Note that here we set the scaling factor for the inverse operation of
1049 // getting the "true" objective value from the scaled one. Hence the
1050 // inverse.
1051 auto* objective_proto = cp_model->mutable_objective();
1052 const int64_t mult = maximize ? -1 : 1;
1053 objective_proto->set_offset(objective_offset * scaling_factor / gcd * mult);
1054 objective_proto->set_scaling_factor(1.0 / scaling_factor * gcd * mult);
1055 for (int i = 0; i < coefficients.size(); ++i) {
1056 const int64_t value =
1057 static_cast<int64_t>(std::round(coefficients[i] * scaling_factor)) /
1058 gcd;
1059 if (value != 0) {
1060 objective_proto->add_vars(var_indices[i]);
1061 objective_proto->add_coeffs(value * mult);
1062 }
1063 }
1064
1065 if (scaled_sum_error == 0.0) {
1066 objective_proto->set_scaling_was_exact(true);
1067 }
1068
1069 return true;
1070 }
1071
ConvertBinaryMPModelProtoToBooleanProblem(const MPModelProto & mp_model,LinearBooleanProblem * problem)1072 bool ConvertBinaryMPModelProtoToBooleanProblem(const MPModelProto& mp_model,
1073 LinearBooleanProblem* problem) {
1074 CHECK(problem != nullptr);
1075 problem->Clear();
1076 problem->set_name(mp_model.name());
1077 const int num_variables = mp_model.variable_size();
1078 problem->set_num_variables(num_variables);
1079
1080 // Test if the variables are binary variables.
1081 // Add constraints for the fixed variables.
1082 for (int var_id(0); var_id < num_variables; ++var_id) {
1083 const MPVariableProto& mp_var = mp_model.variable(var_id);
1084 problem->add_var_names(mp_var.name());
1085
1086 // This will be changed to false as soon as we detect the variable to be
1087 // non-binary. This is done this way so we can display a nice error message
1088 // before aborting the function and returning false.
1089 bool is_binary = mp_var.is_integer();
1090
1091 const Fractional lb = mp_var.lower_bound();
1092 const Fractional ub = mp_var.upper_bound();
1093 if (lb <= -1.0) is_binary = false;
1094 if (ub >= 2.0) is_binary = false;
1095 if (is_binary) {
1096 // 4 cases.
1097 if (lb <= 0.0 && ub >= 1.0) {
1098 // Binary variable. Ok.
1099 } else if (lb <= 1.0 && ub >= 1.0) {
1100 // Fixed variable at 1.
1101 LinearBooleanConstraint* constraint = problem->add_constraints();
1102 constraint->set_lower_bound(1);
1103 constraint->set_upper_bound(1);
1104 constraint->add_literals(var_id + 1);
1105 constraint->add_coefficients(1);
1106 } else if (lb <= 0.0 && ub >= 0.0) {
1107 // Fixed variable at 0.
1108 LinearBooleanConstraint* constraint = problem->add_constraints();
1109 constraint->set_lower_bound(0);
1110 constraint->set_upper_bound(0);
1111 constraint->add_literals(var_id + 1);
1112 constraint->add_coefficients(1);
1113 } else {
1114 // No possible integer value!
1115 is_binary = false;
1116 }
1117 }
1118
1119 // Abort if the variable is not binary.
1120 if (!is_binary) {
1121 LOG(WARNING) << "The variable #" << var_id << " with name "
1122 << mp_var.name() << " is not binary. "
1123 << "lb: " << lb << " ub: " << ub;
1124 return false;
1125 }
1126 }
1127
1128 // Variables needed to scale the double coefficients into int64_t.
1129 const int64_t kInt64Max = std::numeric_limits<int64_t>::max();
1130 double max_relative_error = 0.0;
1131 double max_bound_error = 0.0;
1132 double max_scaling_factor = 0.0;
1133 double relative_error = 0.0;
1134 double scaling_factor = 0.0;
1135 std::vector<double> coefficients;
1136
1137 // Add all constraints.
1138 for (const MPConstraintProto& mp_constraint : mp_model.constraint()) {
1139 LinearBooleanConstraint* constraint = problem->add_constraints();
1140 constraint->set_name(mp_constraint.name());
1141
1142 // First scale the coefficients of the constraints.
1143 coefficients.clear();
1144 const int num_coeffs = mp_constraint.coefficient_size();
1145 for (int i = 0; i < num_coeffs; ++i) {
1146 coefficients.push_back(mp_constraint.coefficient(i));
1147 }
1148 GetBestScalingOfDoublesToInt64(coefficients, kInt64Max, &scaling_factor,
1149 &relative_error);
1150 const int64_t gcd =
1151 ComputeGcdOfRoundedDoubles(coefficients, scaling_factor);
1152 max_relative_error = std::max(relative_error, max_relative_error);
1153 max_scaling_factor = std::max(scaling_factor / gcd, max_scaling_factor);
1154
1155 double bound_error = 0.0;
1156 for (int i = 0; i < num_coeffs; ++i) {
1157 const double scaled_value = mp_constraint.coefficient(i) * scaling_factor;
1158 bound_error += std::abs(round(scaled_value) - scaled_value);
1159 const int64_t value = static_cast<int64_t>(round(scaled_value)) / gcd;
1160 if (value != 0) {
1161 constraint->add_literals(mp_constraint.var_index(i) + 1);
1162 constraint->add_coefficients(value);
1163 }
1164 }
1165 max_bound_error = std::max(max_bound_error, bound_error);
1166
1167 // Add the bounds. Note that we do not pass them to
1168 // GetBestScalingOfDoublesToInt64() because we know that the sum of absolute
1169 // coefficients of the constraint fit on an int64_t. If one of the scaled
1170 // bound overflows, we don't care by how much because in this case the
1171 // constraint is just trivial or unsatisfiable.
1172 const Fractional lb = mp_constraint.lower_bound();
1173 if (lb != -kInfinity) {
1174 if (lb * scaling_factor > static_cast<double>(kInt64Max)) {
1175 LOG(WARNING) << "A constraint is trivially unsatisfiable.";
1176 return false;
1177 }
1178 if (lb * scaling_factor > -static_cast<double>(kInt64Max)) {
1179 // Otherwise, the constraint is not needed.
1180 constraint->set_lower_bound(
1181 static_cast<int64_t>(round(lb * scaling_factor - bound_error)) /
1182 gcd);
1183 }
1184 }
1185 const Fractional ub = mp_constraint.upper_bound();
1186 if (ub != kInfinity) {
1187 if (ub * scaling_factor < -static_cast<double>(kInt64Max)) {
1188 LOG(WARNING) << "A constraint is trivially unsatisfiable.";
1189 return false;
1190 }
1191 if (ub * scaling_factor < static_cast<double>(kInt64Max)) {
1192 // Otherwise, the constraint is not needed.
1193 constraint->set_upper_bound(
1194 static_cast<int64_t>(round(ub * scaling_factor + bound_error)) /
1195 gcd);
1196 }
1197 }
1198 }
1199
1200 // Display the error/scaling without taking into account the objective first.
1201 LOG(INFO) << "Maximum constraint relative error: " << max_relative_error;
1202 LOG(INFO) << "Maximum constraint bound error: " << max_bound_error;
1203 LOG(INFO) << "Maximum constraint scaling factor: " << max_scaling_factor;
1204
1205 // Add the objective.
1206 coefficients.clear();
1207 for (int var_id = 0; var_id < num_variables; ++var_id) {
1208 const MPVariableProto& mp_var = mp_model.variable(var_id);
1209 coefficients.push_back(mp_var.objective_coefficient());
1210 }
1211 GetBestScalingOfDoublesToInt64(coefficients, kInt64Max, &scaling_factor,
1212 &relative_error);
1213 const int64_t gcd = ComputeGcdOfRoundedDoubles(coefficients, scaling_factor);
1214 max_relative_error = std::max(relative_error, max_relative_error);
1215
1216 // Display the objective error/scaling.
1217 LOG(INFO) << "objective relative error: " << relative_error;
1218 LOG(INFO) << "objective scaling factor: " << scaling_factor / gcd;
1219
1220 LinearObjective* objective = problem->mutable_objective();
1221 objective->set_offset(mp_model.objective_offset() * scaling_factor / gcd);
1222
1223 // Note that here we set the scaling factor for the inverse operation of
1224 // getting the "true" objective value from the scaled one. Hence the inverse.
1225 objective->set_scaling_factor(1.0 / scaling_factor * gcd);
1226 for (int var_id = 0; var_id < num_variables; ++var_id) {
1227 const MPVariableProto& mp_var = mp_model.variable(var_id);
1228 const int64_t value =
1229 static_cast<int64_t>(
1230 round(mp_var.objective_coefficient() * scaling_factor)) /
1231 gcd;
1232 if (value != 0) {
1233 objective->add_literals(var_id + 1);
1234 objective->add_coefficients(value);
1235 }
1236 }
1237
1238 // If the problem was a maximization one, we need to modify the objective.
1239 if (mp_model.maximize()) ChangeOptimizationDirection(problem);
1240
1241 // Test the precision of the conversion.
1242 const double kRelativeTolerance = 1e-8;
1243 if (max_relative_error > kRelativeTolerance) {
1244 LOG(WARNING) << "The relative error during double -> int64_t conversion "
1245 << "is too high!";
1246 return false;
1247 }
1248 return true;
1249 }
1250
ConvertBooleanProblemToLinearProgram(const LinearBooleanProblem & problem,glop::LinearProgram * lp)1251 void ConvertBooleanProblemToLinearProgram(const LinearBooleanProblem& problem,
1252 glop::LinearProgram* lp) {
1253 lp->Clear();
1254 for (int i = 0; i < problem.num_variables(); ++i) {
1255 const ColIndex col = lp->CreateNewVariable();
1256 lp->SetVariableType(col, glop::LinearProgram::VariableType::INTEGER);
1257 lp->SetVariableBounds(col, 0.0, 1.0);
1258 }
1259
1260 // Variables name are optional.
1261 if (problem.var_names_size() != 0) {
1262 CHECK_EQ(problem.var_names_size(), problem.num_variables());
1263 for (int i = 0; i < problem.num_variables(); ++i) {
1264 lp->SetVariableName(ColIndex(i), problem.var_names(i));
1265 }
1266 }
1267
1268 for (const LinearBooleanConstraint& constraint : problem.constraints()) {
1269 const RowIndex constraint_index = lp->CreateNewConstraint();
1270 lp->SetConstraintName(constraint_index, constraint.name());
1271 double sum = 0.0;
1272 for (int i = 0; i < constraint.literals_size(); ++i) {
1273 const int literal = constraint.literals(i);
1274 const double coeff = constraint.coefficients(i);
1275 const ColIndex variable_index = ColIndex(abs(literal) - 1);
1276 if (literal < 0) {
1277 sum += coeff;
1278 lp->SetCoefficient(constraint_index, variable_index, -coeff);
1279 } else {
1280 lp->SetCoefficient(constraint_index, variable_index, coeff);
1281 }
1282 }
1283 lp->SetConstraintBounds(
1284 constraint_index,
1285 constraint.has_lower_bound() ? constraint.lower_bound() - sum
1286 : -kInfinity,
1287 constraint.has_upper_bound() ? constraint.upper_bound() - sum
1288 : kInfinity);
1289 }
1290
1291 // Objective.
1292 {
1293 double sum = 0.0;
1294 const LinearObjective& objective = problem.objective();
1295 const double scaling_factor = objective.scaling_factor();
1296 for (int i = 0; i < objective.literals_size(); ++i) {
1297 const int literal = objective.literals(i);
1298 const double coeff =
1299 static_cast<double>(objective.coefficients(i)) * scaling_factor;
1300 const ColIndex variable_index = ColIndex(abs(literal) - 1);
1301 if (literal < 0) {
1302 sum += coeff;
1303 lp->SetObjectiveCoefficient(variable_index, -coeff);
1304 } else {
1305 lp->SetObjectiveCoefficient(variable_index, coeff);
1306 }
1307 }
1308 lp->SetObjectiveOffset((sum + objective.offset()) * scaling_factor);
1309 lp->SetMaximizationProblem(scaling_factor < 0);
1310 }
1311
1312 lp->CleanUp();
1313 }
1314
FixVariablesFromSat(const SatSolver & solver,glop::LinearProgram * lp)1315 int FixVariablesFromSat(const SatSolver& solver, glop::LinearProgram* lp) {
1316 int num_fixed_variables = 0;
1317 const Trail& trail = solver.LiteralTrail();
1318 for (int i = 0; i < trail.Index(); ++i) {
1319 const BooleanVariable var = trail[i].Variable();
1320 const int value = trail[i].IsPositive() ? 1.0 : 0.0;
1321 if (trail.Info(var).level == 0) {
1322 ++num_fixed_variables;
1323 lp->SetVariableBounds(ColIndex(var.value()), value, value);
1324 }
1325 }
1326 return num_fixed_variables;
1327 }
1328
SolveLpAndUseSolutionForSatAssignmentPreference(const glop::LinearProgram & lp,SatSolver * sat_solver,double max_time_in_seconds)1329 bool SolveLpAndUseSolutionForSatAssignmentPreference(
1330 const glop::LinearProgram& lp, SatSolver* sat_solver,
1331 double max_time_in_seconds) {
1332 glop::LPSolver solver;
1333 glop::GlopParameters glop_parameters;
1334 glop_parameters.set_max_time_in_seconds(max_time_in_seconds);
1335 solver.SetParameters(glop_parameters);
1336 const glop::ProblemStatus& status = solver.Solve(lp);
1337 if (status != glop::ProblemStatus::OPTIMAL &&
1338 status != glop::ProblemStatus::IMPRECISE &&
1339 status != glop::ProblemStatus::PRIMAL_FEASIBLE) {
1340 return false;
1341 }
1342 for (ColIndex col(0); col < lp.num_variables(); ++col) {
1343 const Fractional& value = solver.variable_values()[col];
1344 sat_solver->SetAssignmentPreference(
1345 Literal(BooleanVariable(col.value()), round(value) == 1),
1346 1 - std::abs(value - round(value)));
1347 }
1348 return true;
1349 }
1350
SolveLpAndUseIntegerVariableToStartLNS(const glop::LinearProgram & lp,LinearBooleanProblem * problem)1351 bool SolveLpAndUseIntegerVariableToStartLNS(const glop::LinearProgram& lp,
1352 LinearBooleanProblem* problem) {
1353 glop::LPSolver solver;
1354 const glop::ProblemStatus& status = solver.Solve(lp);
1355 if (status != glop::ProblemStatus::OPTIMAL &&
1356 status != glop::ProblemStatus::PRIMAL_FEASIBLE)
1357 return false;
1358 int num_variable_fixed = 0;
1359 for (ColIndex col(0); col < lp.num_variables(); ++col) {
1360 const Fractional tolerance = 1e-5;
1361 const Fractional& value = solver.variable_values()[col];
1362 if (value > 1 - tolerance) {
1363 ++num_variable_fixed;
1364 LinearBooleanConstraint* constraint = problem->add_constraints();
1365 constraint->set_lower_bound(1);
1366 constraint->set_upper_bound(1);
1367 constraint->add_coefficients(1);
1368 constraint->add_literals(col.value() + 1);
1369 } else if (value < tolerance) {
1370 ++num_variable_fixed;
1371 LinearBooleanConstraint* constraint = problem->add_constraints();
1372 constraint->set_lower_bound(0);
1373 constraint->set_upper_bound(0);
1374 constraint->add_coefficients(1);
1375 constraint->add_literals(col.value() + 1);
1376 }
1377 }
1378 LOG(INFO) << "LNS with " << num_variable_fixed << " fixed variables.";
1379 return true;
1380 }
1381
ComputeTrueObjectiveLowerBound(const CpModelProto & model_proto_with_floating_point_objective,const CpObjectiveProto & integer_objective,const int64_t inner_integer_objective_lower_bound)1382 double ComputeTrueObjectiveLowerBound(
1383 const CpModelProto& model_proto_with_floating_point_objective,
1384 const CpObjectiveProto& integer_objective,
1385 const int64_t inner_integer_objective_lower_bound) {
1386 // Create an LP with the correct variable domain.
1387 glop::LinearProgram lp;
1388 const CpModelProto& proto = model_proto_with_floating_point_objective;
1389 for (int i = 0; i < proto.variables().size(); ++i) {
1390 const auto& domain = proto.variables(i).domain();
1391 lp.SetVariableBounds(lp.CreateNewVariable(), static_cast<double>(domain[0]),
1392 static_cast<double>(domain[domain.size() - 1]));
1393 }
1394
1395 // Add the original problem floating point objective.
1396 // This is user given, so we do need to deal with duplicate entries.
1397 const FloatObjectiveProto& float_obj = proto.floating_point_objective();
1398 lp.SetObjectiveOffset(float_obj.offset());
1399 lp.SetMaximizationProblem(float_obj.maximize());
1400 for (int i = 0; i < float_obj.vars().size(); ++i) {
1401 const glop::ColIndex col(float_obj.vars(i));
1402 const double old_value = lp.objective_coefficients()[col];
1403 lp.SetObjectiveCoefficient(col, old_value + float_obj.coeffs(i));
1404 }
1405
1406 // Add a single constraint "integer_objective >= lower_bound".
1407 const glop::RowIndex ct = lp.CreateNewConstraint();
1408 lp.SetConstraintBounds(
1409 ct, static_cast<double>(inner_integer_objective_lower_bound),
1410 std::numeric_limits<double>::infinity());
1411 for (int i = 0; i < integer_objective.vars().size(); ++i) {
1412 lp.SetCoefficient(ct, glop::ColIndex(integer_objective.vars(i)),
1413 static_cast<double>(integer_objective.coeffs(i)));
1414 }
1415
1416 lp.CleanUp();
1417
1418 // This should be fast.
1419 glop::LPSolver solver;
1420 glop::GlopParameters glop_parameters;
1421 glop_parameters.set_max_deterministic_time(10);
1422 glop_parameters.set_change_status_to_imprecise(false);
1423 solver.SetParameters(glop_parameters);
1424 const glop::ProblemStatus& status = solver.Solve(lp);
1425 if (status == glop::ProblemStatus::OPTIMAL) {
1426 return solver.GetObjectiveValue();
1427 }
1428
1429 // Error. Hoperfully this shouldn't happen.
1430 return float_obj.maximize() ? std::numeric_limits<double>::infinity()
1431 : -std::numeric_limits<double>::infinity();
1432 }
1433
1434 } // namespace sat
1435 } // namespace operations_research
1436