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 //
15 
16 #include "ortools/linear_solver/linear_solver.h"
17 
18 #if !defined(_MSC_VER)
19 #include <unistd.h>
20 #endif
21 
22 #include <atomic>
23 #include <cmath>
24 #include <cstddef>
25 #include <cstdint>
26 #include <utility>
27 
28 #include "absl/status/status.h"
29 #include "absl/status/statusor.h"
30 #include "absl/strings/ascii.h"
31 #include "absl/strings/match.h"
32 #include "absl/strings/str_cat.h"
33 #include "absl/strings/str_format.h"
34 #include "absl/strings/str_replace.h"
35 #include "absl/synchronization/mutex.h"
36 #include "absl/synchronization/notification.h"
37 #include "absl/time/time.h"
38 #include "ortools/base/accurate_sum.h"
39 #include "ortools/base/commandlineflags.h"
40 #include "ortools/base/integral_types.h"
41 #include "ortools/base/logging.h"
42 #include "ortools/base/map_util.h"
43 #include "ortools/base/stl_util.h"
44 #include "ortools/base/threadpool.h"
45 #include "ortools/linear_solver/linear_solver.pb.h"
46 #include "ortools/linear_solver/model_exporter.h"
47 #include "ortools/linear_solver/model_validator.h"
48 #include "ortools/port/file.h"
49 #include "ortools/util/fp_utils.h"
50 #include "ortools/util/time_limit.h"
51 
52 ABSL_FLAG(bool, verify_solution, false,
53           "Systematically verify the solution when calling Solve()"
54           ", and change the return value of Solve() to ABNORMAL if"
55           " an error was detected.");
56 ABSL_FLAG(bool, log_verification_errors, true,
57           "If --verify_solution is set: LOG(ERROR) all errors detected"
58           " during the verification of the solution.");
59 ABSL_FLAG(bool, linear_solver_enable_verbose_output, false,
60           "If set, enables verbose output for the solver. Setting this flag"
61           " is the same as calling MPSolver::EnableOutput().");
62 
63 ABSL_FLAG(bool, mpsolver_bypass_model_validation, false,
64           "If set, the user-provided Model won't be verified before Solve()."
65           " Invalid models will typically trigger various error responses"
66           " from the underlying solvers; sometimes crashes.");
67 
68 namespace operations_research {
69 
SolverTypeIsMip(MPModelRequest::SolverType solver_type)70 bool SolverTypeIsMip(MPModelRequest::SolverType solver_type) {
71   switch (solver_type) {
72     case MPModelRequest::GLOP_LINEAR_PROGRAMMING:
73     case MPModelRequest::CLP_LINEAR_PROGRAMMING:
74     case MPModelRequest::GLPK_LINEAR_PROGRAMMING:
75     case MPModelRequest::GUROBI_LINEAR_PROGRAMMING:
76     case MPModelRequest::XPRESS_LINEAR_PROGRAMMING:
77     case MPModelRequest::CPLEX_LINEAR_PROGRAMMING:
78       return false;
79 
80     case MPModelRequest::SCIP_MIXED_INTEGER_PROGRAMMING:
81     case MPModelRequest::GLPK_MIXED_INTEGER_PROGRAMMING:
82     case MPModelRequest::CBC_MIXED_INTEGER_PROGRAMMING:
83     case MPModelRequest::GUROBI_MIXED_INTEGER_PROGRAMMING:
84     case MPModelRequest::KNAPSACK_MIXED_INTEGER_PROGRAMMING:
85     case MPModelRequest::BOP_INTEGER_PROGRAMMING:
86     case MPModelRequest::SAT_INTEGER_PROGRAMMING:
87     case MPModelRequest::XPRESS_MIXED_INTEGER_PROGRAMMING:
88     case MPModelRequest::CPLEX_MIXED_INTEGER_PROGRAMMING:
89       return true;
90   }
91   LOG(DFATAL) << "Invalid SolverType: " << solver_type;
92   return false;
93 }
94 
GetCoefficient(const MPVariable * const var) const95 double MPConstraint::GetCoefficient(const MPVariable* const var) const {
96   DLOG_IF(DFATAL, !interface_->solver_->OwnsVariable(var)) << var;
97   if (var == nullptr) return 0.0;
98   return gtl::FindWithDefault(coefficients_, var, 0.0);
99 }
100 
SetCoefficient(const MPVariable * const var,double coeff)101 void MPConstraint::SetCoefficient(const MPVariable* const var, double coeff) {
102   DLOG_IF(DFATAL, !interface_->solver_->OwnsVariable(var)) << var;
103   if (var == nullptr) return;
104   if (coeff == 0.0) {
105     auto it = coefficients_.find(var);
106     // If setting a coefficient to 0 when this coefficient did not
107     // exist or was already 0, do nothing: skip
108     // interface_->SetCoefficient() and do not store a coefficient in
109     // the map.  Note that if the coefficient being set to 0 did exist
110     // and was not 0, we do have to keep a 0 in the coefficients_ map,
111     // because the extraction of the constraint might rely on it,
112     // depending on the underlying solver.
113     if (it != coefficients_.end() && it->second != 0.0) {
114       const double old_value = it->second;
115       it->second = 0.0;
116       interface_->SetCoefficient(this, var, 0.0, old_value);
117     }
118     return;
119   }
120   auto insertion_result = coefficients_.insert(std::make_pair(var, coeff));
121   const double old_value =
122       insertion_result.second ? 0.0 : insertion_result.first->second;
123   insertion_result.first->second = coeff;
124   interface_->SetCoefficient(this, var, coeff, old_value);
125 }
126 
Clear()127 void MPConstraint::Clear() {
128   interface_->ClearConstraint(this);
129   coefficients_.clear();
130 }
131 
SetBounds(double lb,double ub)132 void MPConstraint::SetBounds(double lb, double ub) {
133   const bool change = lb != lb_ || ub != ub_;
134   lb_ = lb;
135   ub_ = ub;
136   if (change && interface_->constraint_is_extracted(index_)) {
137     interface_->SetConstraintBounds(index_, lb_, ub_);
138   }
139 }
140 
dual_value() const141 double MPConstraint::dual_value() const {
142   if (!interface_->IsContinuous()) {
143     LOG(DFATAL) << "Dual value only available for continuous problems";
144     return 0.0;
145   }
146   if (!interface_->CheckSolutionIsSynchronizedAndExists()) return 0.0;
147   return dual_value_;
148 }
149 
basis_status() const150 MPSolver::BasisStatus MPConstraint::basis_status() const {
151   if (!interface_->IsContinuous()) {
152     LOG(DFATAL) << "Basis status only available for continuous problems";
153     return MPSolver::FREE;
154   }
155   if (!interface_->CheckSolutionIsSynchronizedAndExists()) {
156     return MPSolver::FREE;
157   }
158   // This is done lazily as this method is expected to be rarely used.
159   return interface_->row_status(index_);
160 }
161 
ContainsNewVariables()162 bool MPConstraint::ContainsNewVariables() {
163   const int last_variable_index = interface_->last_variable_index();
164   for (const auto& entry : coefficients_) {
165     const int variable_index = entry.first->index();
166     if (variable_index >= last_variable_index ||
167         !interface_->variable_is_extracted(variable_index)) {
168       return true;
169     }
170   }
171   return false;
172 }
173 
174 // ----- MPObjective -----
175 
GetCoefficient(const MPVariable * const var) const176 double MPObjective::GetCoefficient(const MPVariable* const var) const {
177   DLOG_IF(DFATAL, !interface_->solver_->OwnsVariable(var)) << var;
178   if (var == nullptr) return 0.0;
179   return gtl::FindWithDefault(coefficients_, var, 0.0);
180 }
181 
SetCoefficient(const MPVariable * const var,double coeff)182 void MPObjective::SetCoefficient(const MPVariable* const var, double coeff) {
183   DLOG_IF(DFATAL, !interface_->solver_->OwnsVariable(var)) << var;
184   if (var == nullptr) return;
185   if (coeff == 0.0) {
186     auto it = coefficients_.find(var);
187     // See the discussion on MPConstraint::SetCoefficient() for 0 coefficients,
188     // the same reasoning applies here.
189     if (it == coefficients_.end() || it->second == 0.0) return;
190     it->second = 0.0;
191   } else {
192     coefficients_[var] = coeff;
193   }
194   interface_->SetObjectiveCoefficient(var, coeff);
195 }
196 
SetOffset(double value)197 void MPObjective::SetOffset(double value) {
198   offset_ = value;
199   interface_->SetObjectiveOffset(offset_);
200 }
201 
202 namespace {
CheckLinearExpr(const MPSolver & solver,const LinearExpr & linear_expr)203 void CheckLinearExpr(const MPSolver& solver, const LinearExpr& linear_expr) {
204   for (auto var_value_pair : linear_expr.terms()) {
205     CHECK(solver.OwnsVariable(var_value_pair.first))
206         << "Bad MPVariable* in LinearExpr, did you try adding an integer to an "
207            "MPVariable* directly?";
208   }
209 }
210 }  // namespace
211 
OptimizeLinearExpr(const LinearExpr & linear_expr,bool is_maximization)212 void MPObjective::OptimizeLinearExpr(const LinearExpr& linear_expr,
213                                      bool is_maximization) {
214   CheckLinearExpr(*interface_->solver_, linear_expr);
215   interface_->ClearObjective();
216   coefficients_.clear();
217   SetOffset(linear_expr.offset());
218   for (const auto& kv : linear_expr.terms()) {
219     SetCoefficient(kv.first, kv.second);
220   }
221   SetOptimizationDirection(is_maximization);
222 }
223 
AddLinearExpr(const LinearExpr & linear_expr)224 void MPObjective::AddLinearExpr(const LinearExpr& linear_expr) {
225   CheckLinearExpr(*interface_->solver_, linear_expr);
226   SetOffset(offset_ + linear_expr.offset());
227   for (const auto& kv : linear_expr.terms()) {
228     SetCoefficient(kv.first, GetCoefficient(kv.first) + kv.second);
229   }
230 }
231 
Clear()232 void MPObjective::Clear() {
233   interface_->ClearObjective();
234   coefficients_.clear();
235   offset_ = 0.0;
236   SetMinimization();
237 }
238 
SetOptimizationDirection(bool maximize)239 void MPObjective::SetOptimizationDirection(bool maximize) {
240   // Note(user): The maximize_ bool would more naturally belong to the
241   // MPObjective, but it actually has to be a member of MPSolverInterface,
242   // because some implementations (such as GLPK) need that bool for the
243   // MPSolverInterface constructor, i.e. at a time when the MPObjective is not
244   // constructed yet (MPSolverInterface is always built before MPObjective
245   // when a new MPSolver is constructed).
246   interface_->maximize_ = maximize;
247   interface_->SetOptimizationDirection(maximize);
248 }
249 
maximization() const250 bool MPObjective::maximization() const { return interface_->maximize_; }
251 
minimization() const252 bool MPObjective::minimization() const { return !interface_->maximize_; }
253 
Value() const254 double MPObjective::Value() const {
255   // Note(user): implementation-wise, the objective value belongs more
256   // naturally to the MPSolverInterface, since all of its implementations write
257   // to it directly.
258   return interface_->objective_value();
259 }
260 
BestBound() const261 double MPObjective::BestBound() const {
262   // Note(user): the best objective bound belongs to the interface for the
263   // same reasons as the objective value does.
264   return interface_->best_objective_bound();
265 }
266 
267 // ----- MPVariable -----
268 
solution_value() const269 double MPVariable::solution_value() const {
270   if (!interface_->CheckSolutionIsSynchronizedAndExists()) return 0.0;
271   // If the underlying solver supports integer variables, and this is an integer
272   // variable, we round the solution value (i.e., clients usually expect precise
273   // integer values for integer variables).
274   return (integer_ && interface_->IsMIP()) ? round(solution_value_)
275                                            : solution_value_;
276 }
277 
unrounded_solution_value() const278 double MPVariable::unrounded_solution_value() const {
279   if (!interface_->CheckSolutionIsSynchronizedAndExists()) return 0.0;
280   return solution_value_;
281 }
282 
reduced_cost() const283 double MPVariable::reduced_cost() const {
284   if (!interface_->IsContinuous()) {
285     LOG(DFATAL) << "Reduced cost only available for continuous problems";
286     return 0.0;
287   }
288   if (!interface_->CheckSolutionIsSynchronizedAndExists()) return 0.0;
289   return reduced_cost_;
290 }
291 
basis_status() const292 MPSolver::BasisStatus MPVariable::basis_status() const {
293   if (!interface_->IsContinuous()) {
294     LOG(DFATAL) << "Basis status only available for continuous problems";
295     return MPSolver::FREE;
296   }
297   if (!interface_->CheckSolutionIsSynchronizedAndExists()) {
298     return MPSolver::FREE;
299   }
300   // This is done lazily as this method is expected to be rarely used.
301   return interface_->column_status(index_);
302 }
303 
SetBounds(double lb,double ub)304 void MPVariable::SetBounds(double lb, double ub) {
305   const bool change = lb != lb_ || ub != ub_;
306   lb_ = lb;
307   ub_ = ub;
308   if (change && interface_->variable_is_extracted(index_)) {
309     interface_->SetVariableBounds(index_, lb_, ub_);
310   }
311 }
312 
SetInteger(bool integer)313 void MPVariable::SetInteger(bool integer) {
314   if (integer_ != integer) {
315     integer_ = integer;
316     if (interface_->variable_is_extracted(index_)) {
317       interface_->SetVariableInteger(index_, integer);
318     }
319   }
320 }
321 
SetBranchingPriority(int priority)322 void MPVariable::SetBranchingPriority(int priority) {
323   if (priority == branching_priority_) return;
324   branching_priority_ = priority;
325   interface_->BranchingPriorityChangedForVariable(index_);
326 }
327 
328 // ----- Interface shortcuts -----
329 
IsMIP() const330 bool MPSolver::IsMIP() const { return interface_->IsMIP(); }
331 
SolverVersion() const332 std::string MPSolver::SolverVersion() const {
333   return interface_->SolverVersion();
334 }
335 
underlying_solver()336 void* MPSolver::underlying_solver() { return interface_->underlying_solver(); }
337 
338 // ---- Solver-specific parameters ----
339 
SetNumThreads(int num_threads)340 absl::Status MPSolver::SetNumThreads(int num_threads) {
341   if (num_threads < 1) {
342     return absl::InvalidArgumentError("num_threads must be a positive number.");
343   }
344   const absl::Status status = interface_->SetNumThreads(num_threads);
345   if (status.ok()) {
346     num_threads_ = num_threads;
347   }
348   return status;
349 }
350 
SetSolverSpecificParametersAsString(const std::string & parameters)351 bool MPSolver::SetSolverSpecificParametersAsString(
352     const std::string& parameters) {
353   solver_specific_parameter_string_ = parameters;
354   return interface_->SetSolverSpecificParametersAsString(parameters);
355 }
356 
357 // ----- Solver -----
358 
359 #if defined(USE_CLP) || defined(USE_CBC)
360 extern MPSolverInterface* BuildCLPInterface(MPSolver* const solver);
361 #endif
362 #if defined(USE_CBC)
363 extern MPSolverInterface* BuildCBCInterface(MPSolver* const solver);
364 #endif
365 #if defined(USE_GLPK)
366 extern MPSolverInterface* BuildGLPKInterface(bool mip, MPSolver* const solver);
367 #endif
368 extern MPSolverInterface* BuildBopInterface(MPSolver* const solver);
369 extern MPSolverInterface* BuildGLOPInterface(MPSolver* const solver);
370 extern MPSolverInterface* BuildSatInterface(MPSolver* const solver);
371 #if defined(USE_SCIP)
372 extern MPSolverInterface* BuildSCIPInterface(MPSolver* const solver);
373 #endif
374 extern MPSolverInterface* BuildGurobiInterface(bool mip,
375                                                MPSolver* const solver);
376 #if defined(USE_CPLEX)
377 extern MPSolverInterface* BuildCplexInterface(bool mip, MPSolver* const solver);
378 
379 extern MPSolverInterface* BuildGLOPInterface(MPSolver* const solver);
380 #endif
381 #if defined(USE_XPRESS)
382 extern MPSolverInterface* BuildXpressInterface(bool mip,
383                                                MPSolver* const solver);
384 #endif
385 
386 namespace {
BuildSolverInterface(MPSolver * const solver)387 MPSolverInterface* BuildSolverInterface(MPSolver* const solver) {
388   DCHECK(solver != nullptr);
389   switch (solver->ProblemType()) {
390     case MPSolver::BOP_INTEGER_PROGRAMMING:
391       return BuildBopInterface(solver);
392     case MPSolver::SAT_INTEGER_PROGRAMMING:
393       return BuildSatInterface(solver);
394     case MPSolver::GLOP_LINEAR_PROGRAMMING:
395       return BuildGLOPInterface(solver);
396 #if defined(USE_GLPK)
397     case MPSolver::GLPK_LINEAR_PROGRAMMING:
398       return BuildGLPKInterface(false, solver);
399     case MPSolver::GLPK_MIXED_INTEGER_PROGRAMMING:
400       return BuildGLPKInterface(true, solver);
401 #endif
402 #if defined(USE_CLP) || defined(USE_CBC)
403     case MPSolver::CLP_LINEAR_PROGRAMMING:
404       return BuildCLPInterface(solver);
405 #endif
406 #if defined(USE_CBC)
407     case MPSolver::CBC_MIXED_INTEGER_PROGRAMMING:
408       return BuildCBCInterface(solver);
409 #endif
410 #if defined(USE_SCIP)
411     case MPSolver::SCIP_MIXED_INTEGER_PROGRAMMING:
412       return BuildSCIPInterface(solver);
413 #endif
414     case MPSolver::GUROBI_LINEAR_PROGRAMMING:
415       return BuildGurobiInterface(false, solver);
416     case MPSolver::GUROBI_MIXED_INTEGER_PROGRAMMING:
417       return BuildGurobiInterface(true, solver);
418 #if defined(USE_CPLEX)
419     case MPSolver::CPLEX_LINEAR_PROGRAMMING:
420       return BuildCplexInterface(false, solver);
421     case MPSolver::CPLEX_MIXED_INTEGER_PROGRAMMING:
422       return BuildCplexInterface(true, solver);
423 #endif
424 #if defined(USE_XPRESS)
425     case MPSolver::XPRESS_MIXED_INTEGER_PROGRAMMING:
426       return BuildXpressInterface(true, solver);
427     case MPSolver::XPRESS_LINEAR_PROGRAMMING:
428       return BuildXpressInterface(false, solver);
429 #endif
430     default:
431       // TODO(user): Revert to the best *available* interface.
432       LOG(FATAL) << "Linear solver not recognized.";
433   }
434   return nullptr;
435 }
436 }  // namespace
437 
438 namespace {
NumDigits(int n)439 int NumDigits(int n) {
440 // Number of digits needed to write a non-negative integer in base 10.
441 // Note(user): max(1, log(0) + 1) == max(1, -inf) == 1.
442 #if defined(_MSC_VER)
443   return static_cast<int>(std::max(1.0L, log(1.0L * n) / log(10.0L) + 1.0));
444 #else
445   return static_cast<int>(std::max(1.0, log10(static_cast<double>(n)) + 1.0));
446 #endif
447 }
448 }  // namespace
449 
MPSolver(const std::string & name,OptimizationProblemType problem_type)450 MPSolver::MPSolver(const std::string& name,
451                    OptimizationProblemType problem_type)
452     : name_(name),
453       problem_type_(problem_type),
454       construction_time_(absl::Now()) {
455   interface_.reset(BuildSolverInterface(this));
456   if (absl::GetFlag(FLAGS_linear_solver_enable_verbose_output)) {
457     EnableOutput();
458   }
459   objective_.reset(new MPObjective(interface_.get()));
460 }
461 
~MPSolver()462 MPSolver::~MPSolver() { Clear(); }
463 
464 extern bool GurobiIsCorrectlyInstalled();
465 
466 // static
SupportsProblemType(OptimizationProblemType problem_type)467 bool MPSolver::SupportsProblemType(OptimizationProblemType problem_type) {
468 #ifdef USE_CLP
469   if (problem_type == CLP_LINEAR_PROGRAMMING) return true;
470 #endif
471 #ifdef USE_GLPK
472   if (problem_type == GLPK_LINEAR_PROGRAMMING ||
473       problem_type == GLPK_MIXED_INTEGER_PROGRAMMING) {
474     return true;
475   }
476 #endif
477   if (problem_type == BOP_INTEGER_PROGRAMMING) return true;
478   if (problem_type == SAT_INTEGER_PROGRAMMING) return true;
479   if (problem_type == GLOP_LINEAR_PROGRAMMING) return true;
480   if (problem_type == GUROBI_LINEAR_PROGRAMMING ||
481       problem_type == GUROBI_MIXED_INTEGER_PROGRAMMING) {
482     return GurobiIsCorrectlyInstalled();
483   }
484 #ifdef USE_SCIP
485   if (problem_type == SCIP_MIXED_INTEGER_PROGRAMMING) return true;
486 #endif
487 #ifdef USE_CBC
488   if (problem_type == CBC_MIXED_INTEGER_PROGRAMMING) return true;
489 #endif
490 #ifdef USE_XPRESS
491   if (problem_type == XPRESS_MIXED_INTEGER_PROGRAMMING ||
492       problem_type == XPRESS_LINEAR_PROGRAMMING) {
493     return true;
494   }
495 #endif
496 #ifdef USE_CPLEX
497   if (problem_type == CPLEX_LINEAR_PROGRAMMING ||
498       problem_type == CPLEX_MIXED_INTEGER_PROGRAMMING) {
499     return true;
500   }
501 #endif
502   return false;
503 }
504 
505 // TODO(user): post c++ 14, instead use
506 //   std::pair<MPSolver::OptimizationProblemType, const absl::string_view>
507 // once pair gets a constexpr constructor.
508 namespace {
509 struct NamedOptimizationProblemType {
510   MPSolver::OptimizationProblemType problem_type;
511   absl::string_view name;
512 };
513 }  // namespace
514 
515 #if defined(_MSC_VER)
516 const
517 #else
518 constexpr
519 #endif
520     NamedOptimizationProblemType kOptimizationProblemTypeNames[] = {
521         {MPSolver::GLOP_LINEAR_PROGRAMMING, "glop"},
522         {MPSolver::CLP_LINEAR_PROGRAMMING, "clp"},
523         {MPSolver::GUROBI_LINEAR_PROGRAMMING, "gurobi_lp"},
524         {MPSolver::GLPK_LINEAR_PROGRAMMING, "glpk_lp"},
525         {MPSolver::CPLEX_LINEAR_PROGRAMMING, "cplex_lp"},
526         {MPSolver::XPRESS_LINEAR_PROGRAMMING, "xpress_lp"},
527         {MPSolver::SCIP_MIXED_INTEGER_PROGRAMMING, "scip"},
528         {MPSolver::CBC_MIXED_INTEGER_PROGRAMMING, "cbc"},
529         {MPSolver::SAT_INTEGER_PROGRAMMING, "sat"},
530         {MPSolver::BOP_INTEGER_PROGRAMMING, "bop"},
531         {MPSolver::GUROBI_MIXED_INTEGER_PROGRAMMING, "gurobi"},
532         {MPSolver::GLPK_MIXED_INTEGER_PROGRAMMING, "glpk"},
533         {MPSolver::KNAPSACK_MIXED_INTEGER_PROGRAMMING, "knapsack"},
534         {MPSolver::CPLEX_MIXED_INTEGER_PROGRAMMING, "cplex"},
535         {MPSolver::XPRESS_MIXED_INTEGER_PROGRAMMING, "xpress"},
536 
537 };
538 // static
ParseSolverType(absl::string_view solver_id,MPSolver::OptimizationProblemType * type)539 bool MPSolver::ParseSolverType(absl::string_view solver_id,
540                                MPSolver::OptimizationProblemType* type) {
541   // Normalize the solver id.
542   const std::string id =
543       absl::StrReplaceAll(absl::AsciiStrToUpper(solver_id), {{"-", "_"}});
544 
545   // Support the full enum name
546   MPModelRequest::SolverType solver_type;
547   if (MPModelRequest::SolverType_Parse(id, &solver_type)) {
548     *type = static_cast<MPSolver::OptimizationProblemType>(solver_type);
549     return true;
550   }
551 
552   // Names are stored in lower case.
553   std::string lower_id = absl::AsciiStrToLower(id);
554 
555   // Remove any "_mip" suffix, since they are optional.
556   if (absl::EndsWith(lower_id, "_mip")) {
557     lower_id = lower_id.substr(0, lower_id.size() - 4);
558   }
559 
560   // Rewrite CP-SAT into SAT.
561   if (lower_id == "cp_sat") {
562     lower_id = "sat";
563   }
564 
565   // Reverse lookup in the kOptimizationProblemTypeNames[] array.
566   for (auto& named_solver : kOptimizationProblemTypeNames) {
567     if (named_solver.name == lower_id) {
568       *type = named_solver.problem_type;
569       return true;
570     }
571   }
572 
573   return false;
574 }
575 
ToString(MPSolver::OptimizationProblemType optimization_problem_type)576 const absl::string_view ToString(
577     MPSolver::OptimizationProblemType optimization_problem_type) {
578   for (const auto& named_solver : kOptimizationProblemTypeNames) {
579     if (named_solver.problem_type == optimization_problem_type) {
580       return named_solver.name;
581     }
582   }
583   LOG(FATAL) << "Unrecognized solver type: "
584              << static_cast<int>(optimization_problem_type);
585   return "";
586 }
587 
AbslParseFlag(const absl::string_view text,MPSolver::OptimizationProblemType * solver_type,std::string * error)588 bool AbslParseFlag(const absl::string_view text,
589                    MPSolver::OptimizationProblemType* solver_type,
590                    std::string* error) {
591   DCHECK(solver_type != nullptr);
592   DCHECK(error != nullptr);
593   const bool result = MPSolver::ParseSolverType(text, solver_type);
594   if (!result) {
595     *error = absl::StrCat("Solver type: ", text, " does not exist.");
596   }
597   return result;
598 }
599 
600 /* static */
ParseSolverTypeOrDie(const std::string & solver_id)601 MPSolver::OptimizationProblemType MPSolver::ParseSolverTypeOrDie(
602     const std::string& solver_id) {
603   MPSolver::OptimizationProblemType problem_type;
604   CHECK(MPSolver::ParseSolverType(solver_id, &problem_type)) << solver_id;
605   return problem_type;
606 }
607 
608 /* static */
CreateSolver(const std::string & solver_id)609 MPSolver* MPSolver::CreateSolver(const std::string& solver_id) {
610   MPSolver::OptimizationProblemType problem_type;
611   if (!MPSolver::ParseSolverType(solver_id, &problem_type)) {
612     LOG(WARNING) << "Unrecognized solver type: " << solver_id;
613     return nullptr;
614   }
615   if (!MPSolver::SupportsProblemType(problem_type)) {
616     LOG(WARNING) << "Support for " << solver_id
617                  << " not linked in, or the license was not found.";
618     return nullptr;
619   }
620   MPSolver* solver = new MPSolver("", problem_type);
621   return solver;
622 }
623 
LookupVariableOrNull(const std::string & var_name) const624 MPVariable* MPSolver::LookupVariableOrNull(const std::string& var_name) const {
625   if (!variable_name_to_index_) GenerateVariableNameIndex();
626 
627   absl::flat_hash_map<std::string, int>::const_iterator it =
628       variable_name_to_index_->find(var_name);
629   if (it == variable_name_to_index_->end()) return nullptr;
630   return variables_[it->second];
631 }
632 
LookupConstraintOrNull(const std::string & constraint_name) const633 MPConstraint* MPSolver::LookupConstraintOrNull(
634     const std::string& constraint_name) const {
635   if (!constraint_name_to_index_) GenerateConstraintNameIndex();
636 
637   const auto it = constraint_name_to_index_->find(constraint_name);
638   if (it == constraint_name_to_index_->end()) return nullptr;
639   return constraints_[it->second];
640 }
641 
642 // ----- Methods using protocol buffers -----
643 
LoadModelFromProto(const MPModelProto & input_model,std::string * error_message)644 MPSolverResponseStatus MPSolver::LoadModelFromProto(
645     const MPModelProto& input_model, std::string* error_message) {
646   Clear();
647 
648   // The variable and constraint names are dropped, because we allow
649   // duplicate names in the proto (they're not considered as 'ids'),
650   // unlike the MPSolver C++ API which crashes if there are duplicate names.
651   // Clearing the names makes the MPSolver generate unique names.
652   return LoadModelFromProtoInternal(input_model, /*clear_names=*/true,
653                                     /*check_model_validity=*/true,
654                                     error_message);
655 }
656 
LoadModelFromProtoWithUniqueNamesOrDie(const MPModelProto & input_model,std::string * error_message)657 MPSolverResponseStatus MPSolver::LoadModelFromProtoWithUniqueNamesOrDie(
658     const MPModelProto& input_model, std::string* error_message) {
659   Clear();
660 
661   // Force variable and constraint name indexing (which CHECKs name uniqueness).
662   GenerateVariableNameIndex();
663   GenerateConstraintNameIndex();
664 
665   return LoadModelFromProtoInternal(input_model, /*clear_names=*/false,
666                                     /*check_model_validity=*/true,
667                                     error_message);
668 }
669 
LoadModelFromProtoInternal(const MPModelProto & input_model,bool clear_names,bool check_model_validity,std::string * error_message)670 MPSolverResponseStatus MPSolver::LoadModelFromProtoInternal(
671     const MPModelProto& input_model, bool clear_names,
672     bool check_model_validity, std::string* error_message) {
673   CHECK(error_message != nullptr);
674   if (check_model_validity) {
675     const std::string error = FindErrorInMPModelProto(input_model);
676     if (!error.empty()) {
677       *error_message = error;
678       LOG_IF(INFO, OutputIsEnabled())
679           << "Invalid model given to LoadModelFromProto(): " << error;
680       if (absl::GetFlag(FLAGS_mpsolver_bypass_model_validation)) {
681         LOG_IF(INFO, OutputIsEnabled())
682             << "Ignoring the model error(s) because of"
683             << " --mpsolver_bypass_model_validation.";
684       } else {
685         return absl::StrContains(error, "Infeasible") ? MPSOLVER_INFEASIBLE
686                                                       : MPSOLVER_MODEL_INVALID;
687       }
688     }
689   }
690 
691   if (input_model.has_quadratic_objective()) {
692     *error_message =
693         "Optimizing a quadratic objective is only supported through direct "
694         "proto solves. Please use MPSolver::SolveWithProto, or the solver's "
695         "direct proto solve function.";
696     return MPSOLVER_MODEL_INVALID;
697   }
698 
699   MPObjective* const objective = MutableObjective();
700   // Passing empty names makes the MPSolver generate unique names.
701   const std::string empty;
702   for (int i = 0; i < input_model.variable_size(); ++i) {
703     const MPVariableProto& var_proto = input_model.variable(i);
704     MPVariable* variable =
705         MakeNumVar(var_proto.lower_bound(), var_proto.upper_bound(),
706                    clear_names ? empty : var_proto.name());
707     variable->SetInteger(var_proto.is_integer());
708     if (var_proto.branching_priority() != 0) {
709       variable->SetBranchingPriority(var_proto.branching_priority());
710     }
711     objective->SetCoefficient(variable, var_proto.objective_coefficient());
712   }
713 
714   for (const MPConstraintProto& ct_proto : input_model.constraint()) {
715     if (ct_proto.lower_bound() == -infinity() &&
716         ct_proto.upper_bound() == infinity()) {
717       continue;
718     }
719 
720     MPConstraint* const ct =
721         MakeRowConstraint(ct_proto.lower_bound(), ct_proto.upper_bound(),
722                           clear_names ? empty : ct_proto.name());
723     ct->set_is_lazy(ct_proto.is_lazy());
724     for (int j = 0; j < ct_proto.var_index_size(); ++j) {
725       ct->SetCoefficient(variables_[ct_proto.var_index(j)],
726                          ct_proto.coefficient(j));
727     }
728   }
729 
730   for (const MPGeneralConstraintProto& general_constraint :
731        input_model.general_constraint()) {
732     switch (general_constraint.general_constraint_case()) {
733       case MPGeneralConstraintProto::kIndicatorConstraint: {
734         const auto& proto =
735             general_constraint.indicator_constraint().constraint();
736         if (proto.lower_bound() == -infinity() &&
737             proto.upper_bound() == infinity()) {
738           continue;
739         }
740 
741         const int constraint_index = NumConstraints();
742         MPConstraint* const constraint = new MPConstraint(
743             constraint_index, proto.lower_bound(), proto.upper_bound(),
744             clear_names ? "" : proto.name(), interface_.get());
745         if (constraint_name_to_index_) {
746           gtl::InsertOrDie(&*constraint_name_to_index_, constraint->name(),
747                            constraint_index);
748         }
749         constraints_.push_back(constraint);
750         constraint_is_extracted_.push_back(false);
751 
752         constraint->set_is_lazy(proto.is_lazy());
753         for (int j = 0; j < proto.var_index_size(); ++j) {
754           constraint->SetCoefficient(variables_[proto.var_index(j)],
755                                      proto.coefficient(j));
756         }
757 
758         MPVariable* const variable =
759             variables_[general_constraint.indicator_constraint().var_index()];
760         constraint->indicator_variable_ = variable;
761         constraint->indicator_value_ =
762             general_constraint.indicator_constraint().var_value();
763 
764         if (!interface_->AddIndicatorConstraint(constraint)) {
765           *error_message = "Solver doesn't support indicator constraints";
766           return MPSOLVER_MODEL_INVALID;
767         }
768         break;
769       }
770       default:
771         *error_message = absl::StrFormat(
772             "Optimizing general constraints of type %i is only supported "
773             "through direct proto solves. Please use MPSolver::SolveWithProto, "
774             "or the solver's direct proto solve function.",
775             general_constraint.general_constraint_case());
776         return MPSOLVER_MODEL_INVALID;
777     }
778   }
779 
780   objective->SetOptimizationDirection(input_model.maximize());
781   if (input_model.has_objective_offset()) {
782     objective->SetOffset(input_model.objective_offset());
783   }
784 
785   // Stores any hints about where to start the solve.
786   solution_hint_.clear();
787   for (int i = 0; i < input_model.solution_hint().var_index_size(); ++i) {
788     solution_hint_.push_back(
789         std::make_pair(variables_[input_model.solution_hint().var_index(i)],
790                        input_model.solution_hint().var_value(i)));
791   }
792   return MPSOLVER_MODEL_IS_VALID;
793 }
794 
795 namespace {
ResultStatusToMPSolverResponseStatus(MPSolver::ResultStatus status)796 MPSolverResponseStatus ResultStatusToMPSolverResponseStatus(
797     MPSolver::ResultStatus status) {
798   switch (status) {
799     case MPSolver::OPTIMAL:
800       return MPSOLVER_OPTIMAL;
801     case MPSolver::FEASIBLE:
802       return MPSOLVER_FEASIBLE;
803     case MPSolver::INFEASIBLE:
804       return MPSOLVER_INFEASIBLE;
805     case MPSolver::UNBOUNDED:
806       return MPSOLVER_UNBOUNDED;
807     case MPSolver::ABNORMAL:
808       return MPSOLVER_ABNORMAL;
809     case MPSolver::MODEL_INVALID:
810       return MPSOLVER_MODEL_INVALID;
811     case MPSolver::NOT_SOLVED:
812       return MPSOLVER_NOT_SOLVED;
813   }
814   return MPSOLVER_UNKNOWN_STATUS;
815 }
816 }  // namespace
817 
FillSolutionResponseProto(MPSolutionResponse * response) const818 void MPSolver::FillSolutionResponseProto(MPSolutionResponse* response) const {
819   CHECK(response != nullptr);
820   response->Clear();
821   response->set_status(
822       ResultStatusToMPSolverResponseStatus(interface_->result_status_));
823   if (interface_->result_status_ == MPSolver::OPTIMAL ||
824       interface_->result_status_ == MPSolver::FEASIBLE) {
825     response->set_objective_value(Objective().Value());
826     for (int i = 0; i < variables_.size(); ++i) {
827       response->add_variable_value(variables_[i]->solution_value());
828     }
829 
830     if (interface_->IsMIP()) {
831       response->set_best_objective_bound(interface_->best_objective_bound());
832     } else {
833       // Dual values have no meaning in MIP.
834       for (int j = 0; j < constraints_.size(); ++j) {
835         response->add_dual_value(constraints_[j]->dual_value());
836       }
837       // Reduced cost have no meaning in MIP.
838       for (int i = 0; i < variables_.size(); ++i) {
839         response->add_reduced_cost(variables_[i]->reduced_cost());
840       }
841     }
842   }
843 }
844 
845 namespace {
InCategory(int status,int category)846 bool InCategory(int status, int category) {
847   if (category == MPSOLVER_OPTIMAL) return status == MPSOLVER_OPTIMAL;
848   while (status > category) status >>= 4;
849   return status == category;
850 }
851 
AppendStatusStr(const std::string & msg,MPSolutionResponse * response)852 void AppendStatusStr(const std::string& msg, MPSolutionResponse* response) {
853   response->set_status_str(
854       absl::StrCat(response->status_str(),
855                    (response->status_str().empty() ? "" : "\n"), msg));
856 }
857 }  // namespace
858 
859 // static
SolveWithProto(const MPModelRequest & model_request,MPSolutionResponse * response,std::atomic<bool> * interrupt)860 void MPSolver::SolveWithProto(const MPModelRequest& model_request,
861                               MPSolutionResponse* response,
862                               std::atomic<bool>* interrupt) {
863   CHECK(response != nullptr);
864 
865   if (interrupt != nullptr &&
866       !SolverTypeSupportsInterruption(model_request.solver_type())) {
867     response->set_status(MPSOLVER_INCOMPATIBLE_OPTIONS);
868     response->set_status_str(
869         "Called MPSolver::SolveWithProto with an underlying solver that "
870         "doesn't support interruption.");
871     return;
872   }
873 
874   MPSolver solver(model_request.model().name(),
875                   static_cast<MPSolver::OptimizationProblemType>(
876                       model_request.solver_type()));
877   if (model_request.enable_internal_solver_output()) {
878     solver.EnableOutput();
879   }
880 
881   // If interruption support is not required, we don't need access to the
882   // underlying solver and can solve it directly if the interface supports it.
883   auto optional_response =
884       solver.interface_->DirectlySolveProto(model_request, interrupt);
885   if (optional_response) {
886     *response = std::move(optional_response).value();
887     return;
888   }
889 
890   const absl::optional<LazyMutableCopy<MPModelProto>> optional_model =
891       ExtractValidMPModelOrPopulateResponseStatus(model_request, response);
892   if (!optional_model) {
893     LOG_IF(WARNING, model_request.enable_internal_solver_output())
894         << "Failed to extract a valid model from protocol buffer. Status: "
895         << ProtoEnumToString<MPSolverResponseStatus>(response->status()) << " ("
896         << response->status() << "): " << response->status_str();
897     return;
898   }
899   std::string error_message;
900   response->set_status(solver.LoadModelFromProtoInternal(
901       optional_model->get(), /*clear_names=*/true,
902       /*check_model_validity=*/false, &error_message));
903   // Even though we don't re-check model validity here, there can be some
904   // problems found by LoadModelFromProto, eg. unsupported features.
905   if (response->status() != MPSOLVER_MODEL_IS_VALID) {
906     response->set_status_str(error_message);
907     LOG_IF(WARNING, model_request.enable_internal_solver_output())
908         << "LoadModelFromProtoInternal() failed even though the model was "
909         << "valid! Status: "
910         << ProtoEnumToString<MPSolverResponseStatus>(response->status()) << " ("
911         << response->status() << "); Error: " << error_message;
912     return;
913   }
914   if (model_request.has_solver_time_limit_seconds()) {
915     solver.SetTimeLimit(
916         absl::Seconds(model_request.solver_time_limit_seconds()));
917   }
918   std::string warning_message;
919   if (model_request.has_solver_specific_parameters()) {
920     if (!solver.SetSolverSpecificParametersAsString(
921             model_request.solver_specific_parameters())) {
922       if (model_request.ignore_solver_specific_parameters_failure()) {
923         // We'll add a warning message in status_str after the solve.
924         warning_message =
925             "Warning: the solver specific parameters were not successfully "
926             "applied";
927       } else {
928         response->set_status(MPSOLVER_MODEL_INVALID_SOLVER_PARAMETERS);
929         return;
930       }
931     }
932   }
933 
934   if (interrupt == nullptr) {
935     // If we don't need interruption support, we can save some overhead by
936     // running the solve in the current thread.
937     solver.Solve();
938     solver.FillSolutionResponseProto(response);
939   } else {
940     const absl::Time start_time = absl::Now();
941     absl::Time interrupt_time;
942     bool interrupted_by_user = false;
943     {
944       absl::Notification solve_finished;
945       auto polling_func = [&interrupt, &solve_finished, &solver,
946                            &interrupted_by_user, &interrupt_time,
947                            &model_request]() {
948         constexpr absl::Duration kPollDelay = absl::Microseconds(100);
949         constexpr absl::Duration kMaxInterruptionDelay = absl::Seconds(10);
950 
951         while (!interrupt->load()) {
952           if (solve_finished.HasBeenNotified()) return;
953           absl::SleepFor(kPollDelay);
954         }
955 
956         // If we get here, we received an interruption notification before the
957         // solve finished "naturally".
958         solver.InterruptSolve();
959         interrupt_time = absl::Now();
960         interrupted_by_user = true;
961 
962         // SUBTLE: our call to InterruptSolve() can be ignored by the
963         // underlying solver for several reasons:
964         // 1) The solver thread doesn't poll its 'interrupted' bit often
965         //    enough and takes too long to realize that it should return, or
966         //    its mere return + FillSolutionResponse() takes too long.
967         // 2) The user interrupted the solve so early that Solve() hadn't
968         //    really started yet when we called InterruptSolve().
969         // In case 1), we should just wait a little longer. In case 2), we
970         // should call InterruptSolve() again, maybe several times. To both
971         // accommodate cases where the solver takes really a long time to
972         // react to the interruption, while returning as quickly as possible,
973         // we poll the solve_finished notification with increasing durations
974         // and call InterruptSolve again, each time.
975         for (absl::Duration poll_delay = kPollDelay;
976              absl::Now() <= interrupt_time + kMaxInterruptionDelay;
977              poll_delay *= 2) {
978           if (solve_finished.WaitForNotificationWithTimeout(poll_delay)) {
979             return;
980           } else {
981             solver.InterruptSolve();
982           }
983         }
984 
985         LOG(DFATAL)
986             << "MPSolver::InterruptSolve() seems to be ignored by the "
987                "underlying solver, despite repeated calls over at least "
988             << absl::FormatDuration(kMaxInterruptionDelay)
989             << ". Solver type used: "
990             << MPModelRequest_SolverType_Name(model_request.solver_type());
991 
992         // Note that in opt builds, the polling thread terminates here with an
993         // error message, but we let Solve() finish, ignoring the user
994         // interruption request.
995       };
996 
997       // The choice to do polling rather than solving in the second thread is
998       // not arbitrary, as we want to maintain any custom thread options set by
999       // the user. They shouldn't matter for polling, but for solving we might
1000       // e.g. use a larger stack.
1001       ThreadPool thread_pool("SolverThread", /*num_threads=*/1);
1002       thread_pool.StartWorkers();
1003       thread_pool.Schedule(polling_func);
1004 
1005       // Make sure the interruption notification didn't arrive while waiting to
1006       // be scheduled.
1007       if (!interrupt->load()) {
1008         solver.Solve();
1009         solver.FillSolutionResponseProto(response);
1010       } else {  // *interrupt == true
1011         response->set_status(MPSOLVER_CANCELLED_BY_USER);
1012         response->set_status_str(
1013             "Solve not started, because the user set the atomic<bool> in "
1014             "MPSolver::SolveWithProto() to true before solving could "
1015             "start.");
1016       }
1017       solve_finished.Notify();
1018 
1019       // We block until the thread finishes when thread_pool goes out of scope.
1020     }
1021 
1022     if (interrupted_by_user) {
1023       // Despite the interruption, the solver might still have found a useful
1024       // result. If so, don't overwrite the status.
1025       if (InCategory(response->status(), MPSOLVER_NOT_SOLVED)) {
1026         response->set_status(MPSOLVER_CANCELLED_BY_USER);
1027       }
1028       AppendStatusStr(
1029           absl::StrFormat(
1030               "User interrupted MPSolver::SolveWithProto() by setting the "
1031               "atomic<bool> to true at %s (%s after solving started.)",
1032               absl::FormatTime(interrupt_time),
1033               absl::FormatDuration(interrupt_time - start_time)),
1034           response);
1035     }
1036   }
1037 
1038   if (!warning_message.empty()) {
1039     AppendStatusStr(warning_message, response);
1040   }
1041 }
1042 
ExportModelToProto(MPModelProto * output_model) const1043 void MPSolver::ExportModelToProto(MPModelProto* output_model) const {
1044   DCHECK(output_model != nullptr);
1045   output_model->Clear();
1046   // Name
1047   output_model->set_name(Name());
1048   // Variables
1049   for (int j = 0; j < variables_.size(); ++j) {
1050     const MPVariable* const var = variables_[j];
1051     MPVariableProto* const variable_proto = output_model->add_variable();
1052     // TODO(user): Add option to avoid filling the var name to avoid overly
1053     // large protocol buffers.
1054     variable_proto->set_name(var->name());
1055     variable_proto->set_lower_bound(var->lb());
1056     variable_proto->set_upper_bound(var->ub());
1057     variable_proto->set_is_integer(var->integer());
1058     if (objective_->GetCoefficient(var) != 0.0) {
1059       variable_proto->set_objective_coefficient(
1060           objective_->GetCoefficient(var));
1061     }
1062     if (var->branching_priority() != 0) {
1063       variable_proto->set_branching_priority(var->branching_priority());
1064     }
1065   }
1066 
1067   // Map the variables to their indices. This is needed to output the
1068   // variables in the order they were created, which in turn is needed to have
1069   // repeatable results with ExportModelAsLpFormat and ExportModelAsMpsFormat.
1070   // This step is needed as long as the variable indices are given by the
1071   // underlying solver at the time of model extraction.
1072   // TODO(user): remove this step.
1073   absl::flat_hash_map<const MPVariable*, int> var_to_index;
1074   for (int j = 0; j < variables_.size(); ++j) {
1075     var_to_index[variables_[j]] = j;
1076   }
1077 
1078   // Constraints
1079   for (int i = 0; i < constraints_.size(); ++i) {
1080     MPConstraint* const constraint = constraints_[i];
1081     MPConstraintProto* constraint_proto;
1082     if (constraint->indicator_variable() != nullptr) {
1083       MPGeneralConstraintProto* const general_constraint_proto =
1084           output_model->add_general_constraint();
1085       general_constraint_proto->set_name(constraint->name());
1086       MPIndicatorConstraint* const indicator_constraint_proto =
1087           general_constraint_proto->mutable_indicator_constraint();
1088       indicator_constraint_proto->set_var_index(
1089           constraint->indicator_variable()->index());
1090       indicator_constraint_proto->set_var_value(constraint->indicator_value());
1091       constraint_proto = indicator_constraint_proto->mutable_constraint();
1092     } else {
1093       constraint_proto = output_model->add_constraint();
1094     }
1095     constraint_proto->set_name(constraint->name());
1096     constraint_proto->set_lower_bound(constraint->lb());
1097     constraint_proto->set_upper_bound(constraint->ub());
1098     constraint_proto->set_is_lazy(constraint->is_lazy());
1099     // Vector linear_term will contain pairs (variable index, coeff), that will
1100     // be sorted by variable index.
1101     std::vector<std::pair<int, double>> linear_term;
1102     for (const auto& entry : constraint->coefficients_) {
1103       const MPVariable* const var = entry.first;
1104       const int var_index = gtl::FindWithDefault(var_to_index, var, -1);
1105       DCHECK_NE(-1, var_index);
1106       const double coeff = entry.second;
1107       linear_term.push_back(std::pair<int, double>(var_index, coeff));
1108     }
1109     // The cost of sort is expected to be low as constraints usually have very
1110     // few terms.
1111     std::sort(linear_term.begin(), linear_term.end());
1112     // Now use linear term.
1113     for (const std::pair<int, double>& var_and_coeff : linear_term) {
1114       constraint_proto->add_var_index(var_and_coeff.first);
1115       constraint_proto->add_coefficient(var_and_coeff.second);
1116     }
1117   }
1118 
1119   output_model->set_maximize(Objective().maximization());
1120   output_model->set_objective_offset(Objective().offset());
1121 
1122   if (!solution_hint_.empty()) {
1123     PartialVariableAssignment* const hint =
1124         output_model->mutable_solution_hint();
1125     for (const auto& var_value_pair : solution_hint_) {
1126       hint->add_var_index(var_value_pair.first->index());
1127       hint->add_var_value(var_value_pair.second);
1128     }
1129   }
1130 }
1131 
LoadSolutionFromProto(const MPSolutionResponse & response,double tolerance)1132 absl::Status MPSolver::LoadSolutionFromProto(const MPSolutionResponse& response,
1133                                              double tolerance) {
1134   interface_->result_status_ = static_cast<ResultStatus>(response.status());
1135   if (response.status() != MPSOLVER_OPTIMAL &&
1136       response.status() != MPSOLVER_FEASIBLE) {
1137     return absl::InvalidArgumentError(absl::StrCat(
1138         "Cannot load a solution unless its status is OPTIMAL or FEASIBLE"
1139         " (status was: ",
1140         ProtoEnumToString<MPSolverResponseStatus>(response.status()), ")"));
1141   }
1142   // Before touching the variables, verify that the solution looks legit:
1143   // each variable of the MPSolver must have its value listed exactly once, and
1144   // each listed solution should correspond to a known variable.
1145   if (response.variable_value_size() != variables_.size()) {
1146     return absl::InvalidArgumentError(absl::StrCat(
1147         "Trying to load a solution whose number of variables (",
1148         response.variable_value_size(),
1149         ") does not correspond to the Solver's (", variables_.size(), ")"));
1150   }
1151   interface_->ExtractModel();
1152 
1153   if (tolerance != infinity()) {
1154     // Look further: verify that the variable values are within the bounds.
1155     double largest_error = 0;
1156     int num_vars_out_of_bounds = 0;
1157     int last_offending_var = -1;
1158     for (int i = 0; i < response.variable_value_size(); ++i) {
1159       const double var_value = response.variable_value(i);
1160       MPVariable* var = variables_[i];
1161       // TODO(user): Use parameter when they become available in this class.
1162       const double lb_error = var->lb() - var_value;
1163       const double ub_error = var_value - var->ub();
1164       if (lb_error > tolerance || ub_error > tolerance) {
1165         ++num_vars_out_of_bounds;
1166         largest_error = std::max(largest_error, std::max(lb_error, ub_error));
1167         last_offending_var = i;
1168       }
1169     }
1170     if (num_vars_out_of_bounds > 0) {
1171       return absl::InvalidArgumentError(absl::StrCat(
1172           "Loaded a solution whose variables matched the solver's, but ",
1173           num_vars_out_of_bounds, " of ", variables_.size(),
1174           " variables were out of their bounds, by more than the primal"
1175           " tolerance which is: ",
1176           tolerance, ". Max error: ", largest_error, ", last offender var is #",
1177           last_offending_var, ": '", variables_[last_offending_var]->name(),
1178           "'"));
1179     }
1180   }
1181   for (int i = 0; i < response.variable_value_size(); ++i) {
1182     variables_[i]->set_solution_value(response.variable_value(i));
1183   }
1184   if (response.dual_value_size() > 0) {
1185     if (response.dual_value_size() != constraints_.size()) {
1186       return absl::InvalidArgumentError(absl::StrCat(
1187           "Trying to load a dual solution whose number of entries (",
1188           response.dual_value_size(), ") does not correspond to the Solver's (",
1189           constraints_.size(), ")"));
1190     }
1191     for (int i = 0; i < response.dual_value_size(); ++i) {
1192       constraints_[i]->set_dual_value(response.dual_value(i));
1193     }
1194   }
1195   if (response.reduced_cost_size() > 0) {
1196     if (response.reduced_cost_size() != variables_.size()) {
1197       return absl::InvalidArgumentError(absl::StrCat(
1198           "Trying to load a reduced cost solution whose number of entries (",
1199           response.reduced_cost_size(),
1200           ") does not correspond to the Solver's (", variables_.size(), ")"));
1201     }
1202     for (int i = 0; i < response.reduced_cost_size(); ++i) {
1203       variables_[i]->set_reduced_cost(response.reduced_cost(i));
1204     }
1205   }
1206   // Set the objective value, if is known.
1207   // NOTE(user): We do not verify the objective, even though we could!
1208   if (response.has_objective_value()) {
1209     interface_->objective_value_ = response.objective_value();
1210   }
1211   if (response.has_best_objective_bound()) {
1212     interface_->best_objective_bound_ = response.best_objective_bound();
1213   }
1214   // Mark the status as SOLUTION_SYNCHRONIZED, so that users may inspect the
1215   // solution normally.
1216   interface_->sync_status_ = MPSolverInterface::SOLUTION_SYNCHRONIZED;
1217   return absl::OkStatus();
1218 }
1219 
Clear()1220 void MPSolver::Clear() {
1221   {
1222     absl::MutexLock lock(&global_count_mutex_);
1223     global_num_variables_ += variables_.size();
1224     global_num_constraints_ += constraints_.size();
1225   }
1226   MutableObjective()->Clear();
1227   gtl::STLDeleteElements(&variables_);
1228   gtl::STLDeleteElements(&constraints_);
1229   if (variable_name_to_index_) {
1230     variable_name_to_index_->clear();
1231   }
1232   variable_is_extracted_.clear();
1233   if (constraint_name_to_index_) {
1234     constraint_name_to_index_->clear();
1235   }
1236   constraint_is_extracted_.clear();
1237   interface_->Reset();
1238   solution_hint_.clear();
1239 }
1240 
Reset()1241 void MPSolver::Reset() { interface_->Reset(); }
1242 
InterruptSolve()1243 bool MPSolver::InterruptSolve() { return interface_->InterruptSolve(); }
1244 
SetStartingLpBasis(const std::vector<BasisStatus> & variable_statuses,const std::vector<BasisStatus> & constraint_statuses)1245 void MPSolver::SetStartingLpBasis(
1246     const std::vector<BasisStatus>& variable_statuses,
1247     const std::vector<BasisStatus>& constraint_statuses) {
1248   interface_->SetStartingLpBasis(variable_statuses, constraint_statuses);
1249 }
1250 
MakeVar(double lb,double ub,bool integer,const std::string & name)1251 MPVariable* MPSolver::MakeVar(double lb, double ub, bool integer,
1252                               const std::string& name) {
1253   const int var_index = NumVariables();
1254   MPVariable* v =
1255       new MPVariable(var_index, lb, ub, integer, name, interface_.get());
1256   if (variable_name_to_index_) {
1257     gtl::InsertOrDie(&*variable_name_to_index_, v->name(), var_index);
1258   }
1259   variables_.push_back(v);
1260   variable_is_extracted_.push_back(false);
1261   interface_->AddVariable(v);
1262   return v;
1263 }
1264 
MakeNumVar(double lb,double ub,const std::string & name)1265 MPVariable* MPSolver::MakeNumVar(double lb, double ub,
1266                                  const std::string& name) {
1267   return MakeVar(lb, ub, false, name);
1268 }
1269 
MakeIntVar(double lb,double ub,const std::string & name)1270 MPVariable* MPSolver::MakeIntVar(double lb, double ub,
1271                                  const std::string& name) {
1272   return MakeVar(lb, ub, true, name);
1273 }
1274 
MakeBoolVar(const std::string & name)1275 MPVariable* MPSolver::MakeBoolVar(const std::string& name) {
1276   return MakeVar(0.0, 1.0, true, name);
1277 }
1278 
MakeVarArray(int nb,double lb,double ub,bool integer,const std::string & name,std::vector<MPVariable * > * vars)1279 void MPSolver::MakeVarArray(int nb, double lb, double ub, bool integer,
1280                             const std::string& name,
1281                             std::vector<MPVariable*>* vars) {
1282   DCHECK_GE(nb, 0);
1283   if (nb <= 0) return;
1284   const int num_digits = NumDigits(nb);
1285   for (int i = 0; i < nb; ++i) {
1286     if (name.empty()) {
1287       vars->push_back(MakeVar(lb, ub, integer, name));
1288     } else {
1289       std::string vname =
1290           absl::StrFormat("%s%0*d", name.c_str(), num_digits, i);
1291       vars->push_back(MakeVar(lb, ub, integer, vname));
1292     }
1293   }
1294 }
1295 
MakeNumVarArray(int nb,double lb,double ub,const std::string & name,std::vector<MPVariable * > * vars)1296 void MPSolver::MakeNumVarArray(int nb, double lb, double ub,
1297                                const std::string& name,
1298                                std::vector<MPVariable*>* vars) {
1299   MakeVarArray(nb, lb, ub, false, name, vars);
1300 }
1301 
MakeIntVarArray(int nb,double lb,double ub,const std::string & name,std::vector<MPVariable * > * vars)1302 void MPSolver::MakeIntVarArray(int nb, double lb, double ub,
1303                                const std::string& name,
1304                                std::vector<MPVariable*>* vars) {
1305   MakeVarArray(nb, lb, ub, true, name, vars);
1306 }
1307 
MakeBoolVarArray(int nb,const std::string & name,std::vector<MPVariable * > * vars)1308 void MPSolver::MakeBoolVarArray(int nb, const std::string& name,
1309                                 std::vector<MPVariable*>* vars) {
1310   MakeVarArray(nb, 0.0, 1.0, true, name, vars);
1311 }
1312 
MakeRowConstraint(double lb,double ub)1313 MPConstraint* MPSolver::MakeRowConstraint(double lb, double ub) {
1314   return MakeRowConstraint(lb, ub, "");
1315 }
1316 
MakeRowConstraint()1317 MPConstraint* MPSolver::MakeRowConstraint() {
1318   return MakeRowConstraint(-infinity(), infinity(), "");
1319 }
1320 
MakeRowConstraint(double lb,double ub,const std::string & name)1321 MPConstraint* MPSolver::MakeRowConstraint(double lb, double ub,
1322                                           const std::string& name) {
1323   const int constraint_index = NumConstraints();
1324   MPConstraint* const constraint =
1325       new MPConstraint(constraint_index, lb, ub, name, interface_.get());
1326   if (constraint_name_to_index_) {
1327     gtl::InsertOrDie(&*constraint_name_to_index_, constraint->name(),
1328                      constraint_index);
1329   }
1330   constraints_.push_back(constraint);
1331   constraint_is_extracted_.push_back(false);
1332   interface_->AddRowConstraint(constraint);
1333   return constraint;
1334 }
1335 
MakeRowConstraint(const std::string & name)1336 MPConstraint* MPSolver::MakeRowConstraint(const std::string& name) {
1337   return MakeRowConstraint(-infinity(), infinity(), name);
1338 }
1339 
MakeRowConstraint(const LinearRange & range)1340 MPConstraint* MPSolver::MakeRowConstraint(const LinearRange& range) {
1341   return MakeRowConstraint(range, "");
1342 }
1343 
MakeRowConstraint(const LinearRange & range,const std::string & name)1344 MPConstraint* MPSolver::MakeRowConstraint(const LinearRange& range,
1345                                           const std::string& name) {
1346   CheckLinearExpr(*this, range.linear_expr());
1347   MPConstraint* constraint =
1348       MakeRowConstraint(range.lower_bound(), range.upper_bound(), name);
1349   for (const auto& kv : range.linear_expr().terms()) {
1350     constraint->SetCoefficient(kv.first, kv.second);
1351   }
1352   return constraint;
1353 }
1354 
ComputeMaxConstraintSize(int min_constraint_index,int max_constraint_index) const1355 int MPSolver::ComputeMaxConstraintSize(int min_constraint_index,
1356                                        int max_constraint_index) const {
1357   int max_constraint_size = 0;
1358   DCHECK_GE(min_constraint_index, 0);
1359   DCHECK_LE(max_constraint_index, constraints_.size());
1360   for (int i = min_constraint_index; i < max_constraint_index; ++i) {
1361     MPConstraint* const ct = constraints_[i];
1362     if (ct->coefficients_.size() > max_constraint_size) {
1363       max_constraint_size = ct->coefficients_.size();
1364     }
1365   }
1366   return max_constraint_size;
1367 }
1368 
HasInfeasibleConstraints() const1369 bool MPSolver::HasInfeasibleConstraints() const {
1370   bool hasInfeasibleConstraints = false;
1371   for (int i = 0; i < constraints_.size(); ++i) {
1372     if (constraints_[i]->lb() > constraints_[i]->ub()) {
1373       LOG(WARNING) << "Constraint " << constraints_[i]->name() << " (" << i
1374                    << ") has contradictory bounds:"
1375                    << " lower bound = " << constraints_[i]->lb()
1376                    << " upper bound = " << constraints_[i]->ub();
1377       hasInfeasibleConstraints = true;
1378     }
1379   }
1380   return hasInfeasibleConstraints;
1381 }
1382 
HasIntegerVariables() const1383 bool MPSolver::HasIntegerVariables() const {
1384   for (const MPVariable* const variable : variables_) {
1385     if (variable->integer()) return true;
1386   }
1387   return false;
1388 }
1389 
Solve()1390 MPSolver::ResultStatus MPSolver::Solve() {
1391   MPSolverParameters default_param;
1392   return Solve(default_param);
1393 }
1394 
Solve(const MPSolverParameters & param)1395 MPSolver::ResultStatus MPSolver::Solve(const MPSolverParameters& param) {
1396   // Special case for infeasible constraints so that all solvers have
1397   // the same behavior.
1398   // TODO(user): replace this by model extraction to proto + proto validation
1399   // (the proto has very low overhead compared to the wrapper, both in
1400   // performance and memory, so it's ok).
1401   if (HasInfeasibleConstraints()) {
1402     interface_->result_status_ = MPSolver::INFEASIBLE;
1403     return interface_->result_status_;
1404   }
1405 
1406   MPSolver::ResultStatus status = interface_->Solve(param);
1407   if (absl::GetFlag(FLAGS_verify_solution)) {
1408     if (status != MPSolver::OPTIMAL && status != MPSolver::FEASIBLE) {
1409       VLOG(1) << "--verify_solution enabled, but the solver did not find a"
1410               << " solution: skipping the verification.";
1411     } else if (!VerifySolution(
1412                    param.GetDoubleParam(MPSolverParameters::PRIMAL_TOLERANCE),
1413                    absl::GetFlag(FLAGS_log_verification_errors))) {
1414       status = MPSolver::ABNORMAL;
1415       interface_->result_status_ = status;
1416     }
1417   }
1418   DCHECK_EQ(interface_->result_status_, status);
1419   return status;
1420 }
1421 
Write(const std::string & file_name)1422 void MPSolver::Write(const std::string& file_name) {
1423   interface_->Write(file_name);
1424 }
1425 
1426 namespace {
PrettyPrintVar(const MPVariable & var)1427 std::string PrettyPrintVar(const MPVariable& var) {
1428   const std::string prefix = "Variable '" + var.name() + "': domain = ";
1429   if (var.lb() >= MPSolver::infinity() || var.ub() <= -MPSolver::infinity() ||
1430       var.lb() > var.ub()) {
1431     return prefix + "∅";  // Empty set.
1432   }
1433   // Special case: integer variable with at most two possible values
1434   // (and potentially none).
1435   if (var.integer() && var.ub() - var.lb() <= 1) {
1436     const int64_t lb = static_cast<int64_t>(ceil(var.lb()));
1437     const int64_t ub = static_cast<int64_t>(floor(var.ub()));
1438     if (lb > ub) {
1439       return prefix + "∅";
1440     } else if (lb == ub) {
1441       return absl::StrFormat("%s{ %d }", prefix.c_str(), lb);
1442     } else {
1443       return absl::StrFormat("%s{ %d, %d }", prefix.c_str(), lb, ub);
1444     }
1445   }
1446   // Special case: single (non-infinite) real value.
1447   if (var.lb() == var.ub()) {
1448     return absl::StrFormat("%s{ %f }", prefix.c_str(), var.lb());
1449   }
1450   return prefix + (var.integer() ? "Integer" : "Real") + " in " +
1451          (var.lb() <= -MPSolver::infinity()
1452               ? std::string("]-∞")
1453               : absl::StrFormat("[%f", var.lb())) +
1454          ", " +
1455          (var.ub() >= MPSolver::infinity() ? std::string("+∞[")
1456                                            : absl::StrFormat("%f]", var.ub()));
1457 }
1458 
PrettyPrintConstraint(const MPConstraint & constraint)1459 std::string PrettyPrintConstraint(const MPConstraint& constraint) {
1460   std::string prefix = "Constraint '" + constraint.name() + "': ";
1461   if (constraint.lb() >= MPSolver::infinity() ||
1462       constraint.ub() <= -MPSolver::infinity() ||
1463       constraint.lb() > constraint.ub()) {
1464     return prefix + "ALWAYS FALSE";
1465   }
1466   if (constraint.lb() <= -MPSolver::infinity() &&
1467       constraint.ub() >= MPSolver::infinity()) {
1468     return prefix + "ALWAYS TRUE";
1469   }
1470   prefix += "<linear expr>";
1471   // Equality.
1472   if (constraint.lb() == constraint.ub()) {
1473     return absl::StrFormat("%s = %f", prefix.c_str(), constraint.lb());
1474   }
1475   // Inequalities.
1476   if (constraint.lb() <= -MPSolver::infinity()) {
1477     return absl::StrFormat("%s ≤ %f", prefix.c_str(), constraint.ub());
1478   }
1479   if (constraint.ub() >= MPSolver::infinity()) {
1480     return absl::StrFormat("%s ≥ %f", prefix.c_str(), constraint.lb());
1481   }
1482   return absl::StrFormat("%s ∈ [%f, %f]", prefix.c_str(), constraint.lb(),
1483                          constraint.ub());
1484 }
1485 }  // namespace
1486 
ClampSolutionWithinBounds()1487 absl::Status MPSolver::ClampSolutionWithinBounds() {
1488   interface_->ExtractModel();
1489   for (MPVariable* const variable : variables_) {
1490     const double value = variable->solution_value();
1491     if (std::isnan(value)) {
1492       return absl::InvalidArgumentError(
1493           absl::StrCat("NaN value for ", PrettyPrintVar(*variable)));
1494     }
1495     if (value < variable->lb()) {
1496       variable->set_solution_value(variable->lb());
1497     } else if (value > variable->ub()) {
1498       variable->set_solution_value(variable->ub());
1499     }
1500   }
1501   interface_->sync_status_ = MPSolverInterface::SOLUTION_SYNCHRONIZED;
1502   return absl::OkStatus();
1503 }
1504 
ComputeConstraintActivities() const1505 std::vector<double> MPSolver::ComputeConstraintActivities() const {
1506   // TODO(user): test this failure case.
1507   if (!interface_->CheckSolutionIsSynchronizedAndExists()) return {};
1508   std::vector<double> activities(constraints_.size(), 0.0);
1509   for (int i = 0; i < constraints_.size(); ++i) {
1510     const MPConstraint& constraint = *constraints_[i];
1511     AccurateSum<double> sum;
1512     for (const auto& entry : constraint.coefficients_) {
1513       sum.Add(entry.first->solution_value() * entry.second);
1514     }
1515     activities[i] = sum.Value();
1516   }
1517   return activities;
1518 }
1519 
1520 // TODO(user): split.
VerifySolution(double tolerance,bool log_errors) const1521 bool MPSolver::VerifySolution(double tolerance, bool log_errors) const {
1522   double max_observed_error = 0;
1523   if (tolerance < 0) tolerance = infinity();
1524   int num_errors = 0;
1525 
1526   // Verify variables.
1527   for (int i = 0; i < variables_.size(); ++i) {
1528     const MPVariable& var = *variables_[i];
1529     const double value = var.solution_value();
1530     // Check for NaN.
1531     if (std::isnan(value)) {
1532       ++num_errors;
1533       max_observed_error = infinity();
1534       LOG_IF(ERROR, log_errors) << "NaN value for " << PrettyPrintVar(var);
1535       continue;
1536     }
1537     // Check lower bound.
1538     if (var.lb() != -infinity()) {
1539       if (value < var.lb() - tolerance) {
1540         ++num_errors;
1541         max_observed_error = std::max(max_observed_error, var.lb() - value);
1542         LOG_IF(ERROR, log_errors)
1543             << "Value " << value << " too low for " << PrettyPrintVar(var);
1544       }
1545     }
1546     // Check upper bound.
1547     if (var.ub() != infinity()) {
1548       if (value > var.ub() + tolerance) {
1549         ++num_errors;
1550         max_observed_error = std::max(max_observed_error, value - var.ub());
1551         LOG_IF(ERROR, log_errors)
1552             << "Value " << value << " too high for " << PrettyPrintVar(var);
1553       }
1554     }
1555     // Check integrality.
1556     if (IsMIP() && var.integer()) {
1557       if (fabs(value - round(value)) > tolerance) {
1558         ++num_errors;
1559         max_observed_error =
1560             std::max(max_observed_error, fabs(value - round(value)));
1561         LOG_IF(ERROR, log_errors)
1562             << "Non-integer value " << value << " for " << PrettyPrintVar(var);
1563       }
1564     }
1565   }
1566   if (!IsMIP() && HasIntegerVariables()) {
1567     LOG_IF(INFO, log_errors) << "Skipped variable integrality check, because "
1568                              << "a continuous relaxation of the model was "
1569                              << "solved (i.e., the selected solver does not "
1570                              << "support integer variables).";
1571   }
1572 
1573   // Verify constraints.
1574   const std::vector<double> activities = ComputeConstraintActivities();
1575   for (int i = 0; i < constraints_.size(); ++i) {
1576     const MPConstraint& constraint = *constraints_[i];
1577     const double activity = activities[i];
1578     // Re-compute the activity with a inaccurate summing algorithm.
1579     double inaccurate_activity = 0.0;
1580     for (const auto& entry : constraint.coefficients_) {
1581       inaccurate_activity += entry.first->solution_value() * entry.second;
1582     }
1583     // Catch NaNs.
1584     if (std::isnan(activity) || std::isnan(inaccurate_activity)) {
1585       ++num_errors;
1586       max_observed_error = infinity();
1587       LOG_IF(ERROR, log_errors)
1588           << "NaN value for " << PrettyPrintConstraint(constraint);
1589       continue;
1590     }
1591     // Check bounds.
1592     if (constraint.indicator_variable() == nullptr ||
1593         std::round(constraint.indicator_variable()->solution_value()) ==
1594             constraint.indicator_value()) {
1595       if (constraint.lb() != -infinity()) {
1596         if (activity < constraint.lb() - tolerance) {
1597           ++num_errors;
1598           max_observed_error =
1599               std::max(max_observed_error, constraint.lb() - activity);
1600           LOG_IF(ERROR, log_errors)
1601               << "Activity " << activity << " too low for "
1602               << PrettyPrintConstraint(constraint);
1603         } else if (inaccurate_activity < constraint.lb() - tolerance) {
1604           LOG_IF(WARNING, log_errors)
1605               << "Activity " << activity << ", computed with the (inaccurate)"
1606               << " standard sum of its terms, is too low for "
1607               << PrettyPrintConstraint(constraint);
1608         }
1609       }
1610       if (constraint.ub() != infinity()) {
1611         if (activity > constraint.ub() + tolerance) {
1612           ++num_errors;
1613           max_observed_error =
1614               std::max(max_observed_error, activity - constraint.ub());
1615           LOG_IF(ERROR, log_errors)
1616               << "Activity " << activity << " too high for "
1617               << PrettyPrintConstraint(constraint);
1618         } else if (inaccurate_activity > constraint.ub() + tolerance) {
1619           LOG_IF(WARNING, log_errors)
1620               << "Activity " << activity << ", computed with the (inaccurate)"
1621               << " standard sum of its terms, is too high for "
1622               << PrettyPrintConstraint(constraint);
1623         }
1624       }
1625     }
1626   }
1627 
1628   // Verify that the objective value wasn't reported incorrectly.
1629   const MPObjective& objective = Objective();
1630   AccurateSum<double> objective_sum;
1631   objective_sum.Add(objective.offset());
1632   double inaccurate_objective_value = objective.offset();
1633   for (const auto& entry : objective.coefficients_) {
1634     const double term = entry.first->solution_value() * entry.second;
1635     objective_sum.Add(term);
1636     inaccurate_objective_value += term;
1637   }
1638   const double actual_objective_value = objective_sum.Value();
1639   if (!AreWithinAbsoluteOrRelativeTolerances(
1640           objective.Value(), actual_objective_value, tolerance, tolerance)) {
1641     ++num_errors;
1642     max_observed_error = std::max(
1643         max_observed_error, fabs(actual_objective_value - objective.Value()));
1644     LOG_IF(ERROR, log_errors)
1645         << "Objective value " << objective.Value() << " isn't accurate"
1646         << ", it should be " << actual_objective_value
1647         << " (delta=" << actual_objective_value - objective.Value() << ").";
1648   } else if (!AreWithinAbsoluteOrRelativeTolerances(objective.Value(),
1649                                                     inaccurate_objective_value,
1650                                                     tolerance, tolerance)) {
1651     LOG_IF(WARNING, log_errors)
1652         << "Objective value " << objective.Value() << " doesn't correspond"
1653         << " to the value computed with the standard (and therefore inaccurate)"
1654         << " sum of its terms.";
1655   }
1656   if (num_errors > 0) {
1657     LOG_IF(ERROR, log_errors)
1658         << "There were " << num_errors << " errors above the tolerance ("
1659         << tolerance << "), the largest was " << max_observed_error;
1660     return false;
1661   }
1662   return true;
1663 }
1664 
OutputIsEnabled() const1665 bool MPSolver::OutputIsEnabled() const { return !interface_->quiet(); }
1666 
EnableOutput()1667 void MPSolver::EnableOutput() { interface_->set_quiet(false); }
1668 
SuppressOutput()1669 void MPSolver::SuppressOutput() { interface_->set_quiet(true); }
1670 
iterations() const1671 int64_t MPSolver::iterations() const { return interface_->iterations(); }
1672 
nodes() const1673 int64_t MPSolver::nodes() const { return interface_->nodes(); }
1674 
ComputeExactConditionNumber() const1675 double MPSolver::ComputeExactConditionNumber() const {
1676   return interface_->ComputeExactConditionNumber();
1677 }
1678 
OwnsVariable(const MPVariable * var) const1679 bool MPSolver::OwnsVariable(const MPVariable* var) const {
1680   if (var == nullptr) return false;
1681   if (var->index() >= 0 && var->index() < variables_.size()) {
1682     // Then, verify that the variable with this index has the same address.
1683     return variables_[var->index()] == var;
1684   }
1685   return false;
1686 }
1687 
ExportModelAsLpFormat(bool obfuscate,std::string * model_str) const1688 bool MPSolver::ExportModelAsLpFormat(bool obfuscate,
1689                                      std::string* model_str) const {
1690   MPModelProto proto;
1691   ExportModelToProto(&proto);
1692   MPModelExportOptions options;
1693   options.obfuscate = obfuscate;
1694   const auto status_or =
1695       operations_research::ExportModelAsLpFormat(proto, options);
1696   *model_str = status_or.value_or("");
1697   return status_or.ok();
1698 }
1699 
ExportModelAsMpsFormat(bool fixed_format,bool obfuscate,std::string * model_str) const1700 bool MPSolver::ExportModelAsMpsFormat(bool fixed_format, bool obfuscate,
1701                                       std::string* model_str) const {
1702   MPModelProto proto;
1703   ExportModelToProto(&proto);
1704   MPModelExportOptions options;
1705   options.obfuscate = obfuscate;
1706   const auto status_or =
1707       operations_research::ExportModelAsMpsFormat(proto, options);
1708   *model_str = status_or.value_or("");
1709   return status_or.ok();
1710 }
1711 
SetHint(std::vector<std::pair<const MPVariable *,double>> hint)1712 void MPSolver::SetHint(std::vector<std::pair<const MPVariable*, double>> hint) {
1713   for (const auto& var_value_pair : hint) {
1714     CHECK(OwnsVariable(var_value_pair.first))
1715         << "hint variable does not belong to this solver";
1716   }
1717   solution_hint_ = std::move(hint);
1718 }
1719 
GenerateVariableNameIndex() const1720 void MPSolver::GenerateVariableNameIndex() const {
1721   if (variable_name_to_index_) return;
1722   variable_name_to_index_ = absl::flat_hash_map<std::string, int>();
1723   for (const MPVariable* const var : variables_) {
1724     gtl::InsertOrDie(&*variable_name_to_index_, var->name(), var->index());
1725   }
1726 }
1727 
GenerateConstraintNameIndex() const1728 void MPSolver::GenerateConstraintNameIndex() const {
1729   if (constraint_name_to_index_) return;
1730   constraint_name_to_index_ = absl::flat_hash_map<std::string, int>();
1731   for (const MPConstraint* const cst : constraints_) {
1732     gtl::InsertOrDie(&*constraint_name_to_index_, cst->name(), cst->index());
1733   }
1734 }
1735 
NextSolution()1736 bool MPSolver::NextSolution() { return interface_->NextSolution(); }
1737 
SetCallback(MPCallback * mp_callback)1738 void MPSolver::SetCallback(MPCallback* mp_callback) {
1739   interface_->SetCallback(mp_callback);
1740 }
1741 
SupportsCallbacks() const1742 bool MPSolver::SupportsCallbacks() const {
1743   return interface_->SupportsCallbacks();
1744 }
1745 
1746 // Global counters.
1747 absl::Mutex MPSolver::global_count_mutex_(absl::kConstInit);
1748 int64_t MPSolver::global_num_variables_ = 0;
1749 int64_t MPSolver::global_num_constraints_ = 0;
1750 
1751 // static
global_num_variables()1752 int64_t MPSolver::global_num_variables() {
1753   // Why not ReaderMutexLock? See go/totw/197#when-are-shared-locks-useful.
1754   absl::MutexLock lock(&global_count_mutex_);
1755   return global_num_variables_;
1756 }
1757 
1758 // static
global_num_constraints()1759 int64_t MPSolver::global_num_constraints() {
1760   // Why not ReaderMutexLock? See go/totw/197#when-are-shared-locks-useful.
1761   absl::MutexLock lock(&global_count_mutex_);
1762   return global_num_constraints_;
1763 }
1764 
MPSolverResponseStatusIsRpcError(MPSolverResponseStatus status)1765 bool MPSolverResponseStatusIsRpcError(MPSolverResponseStatus status) {
1766   switch (status) {
1767     // Cases that don't yield an RPC error when they happen on the server.
1768     case MPSOLVER_OPTIMAL:
1769     case MPSOLVER_FEASIBLE:
1770     case MPSOLVER_INFEASIBLE:
1771     case MPSOLVER_NOT_SOLVED:
1772     case MPSOLVER_UNBOUNDED:
1773     case MPSOLVER_ABNORMAL:
1774     case MPSOLVER_UNKNOWN_STATUS:
1775       return false;
1776     // Cases that should never happen with the linear solver server. We prefer
1777     // to consider those as "not RPC errors".
1778     case MPSOLVER_MODEL_IS_VALID:
1779     case MPSOLVER_CANCELLED_BY_USER:
1780       return false;
1781     // Cases that yield an RPC error when they happen on the server.
1782     case MPSOLVER_MODEL_INVALID:
1783     case MPSOLVER_MODEL_INVALID_SOLUTION_HINT:
1784     case MPSOLVER_MODEL_INVALID_SOLVER_PARAMETERS:
1785     case MPSOLVER_SOLVER_TYPE_UNAVAILABLE:
1786     case MPSOLVER_INCOMPATIBLE_OPTIONS:
1787       return true;
1788   }
1789   LOG(DFATAL)
1790       << "MPSolverResponseStatusIsRpcError() called with invalid status "
1791       << "(value: " << status << ")";
1792   return false;
1793 }
1794 
1795 // ---------- MPSolverInterface ----------
1796 
1797 const int MPSolverInterface::kDummyVariableIndex = 0;
1798 
1799 // TODO(user): Initialize objective value and bound to +/- inf (depending on
1800 // optimization direction).
MPSolverInterface(MPSolver * const solver)1801 MPSolverInterface::MPSolverInterface(MPSolver* const solver)
1802     : solver_(solver),
1803       sync_status_(MODEL_SYNCHRONIZED),
1804       result_status_(MPSolver::NOT_SOLVED),
1805       maximize_(false),
1806       last_constraint_index_(0),
1807       last_variable_index_(0),
1808       objective_value_(0.0),
1809       best_objective_bound_(0.0),
1810       quiet_(true) {}
1811 
~MPSolverInterface()1812 MPSolverInterface::~MPSolverInterface() {}
1813 
Write(const std::string & filename)1814 void MPSolverInterface::Write(const std::string& filename) {
1815   LOG(WARNING) << "Writing model not implemented in this solver interface.";
1816 }
1817 
ExtractModel()1818 void MPSolverInterface::ExtractModel() {
1819   switch (sync_status_) {
1820     case MUST_RELOAD: {
1821       ExtractNewVariables();
1822       ExtractNewConstraints();
1823       ExtractObjective();
1824 
1825       last_constraint_index_ = solver_->constraints_.size();
1826       last_variable_index_ = solver_->variables_.size();
1827       sync_status_ = MODEL_SYNCHRONIZED;
1828       break;
1829     }
1830     case MODEL_SYNCHRONIZED: {
1831       // Everything has already been extracted.
1832       DCHECK_EQ(last_constraint_index_, solver_->constraints_.size());
1833       DCHECK_EQ(last_variable_index_, solver_->variables_.size());
1834       break;
1835     }
1836     case SOLUTION_SYNCHRONIZED: {
1837       // Nothing has changed since last solve.
1838       DCHECK_EQ(last_constraint_index_, solver_->constraints_.size());
1839       DCHECK_EQ(last_variable_index_, solver_->variables_.size());
1840       break;
1841     }
1842   }
1843 }
1844 
1845 // TODO(user): remove this method.
ResetExtractionInformation()1846 void MPSolverInterface::ResetExtractionInformation() {
1847   sync_status_ = MUST_RELOAD;
1848   last_constraint_index_ = 0;
1849   last_variable_index_ = 0;
1850   solver_->variable_is_extracted_.assign(solver_->variables_.size(), false);
1851   solver_->constraint_is_extracted_.assign(solver_->constraints_.size(), false);
1852 }
1853 
CheckSolutionIsSynchronized() const1854 bool MPSolverInterface::CheckSolutionIsSynchronized() const {
1855   if (sync_status_ != SOLUTION_SYNCHRONIZED) {
1856     LOG(DFATAL)
1857         << "The model has been changed since the solution was last computed."
1858         << " MPSolverInterface::sync_status_ = " << sync_status_;
1859     return false;
1860   }
1861   return true;
1862 }
1863 
1864 // Default version that can be overwritten by a solver-specific
1865 // version to accommodate for the quirks of each solver.
CheckSolutionExists() const1866 bool MPSolverInterface::CheckSolutionExists() const {
1867   if (result_status_ != MPSolver::OPTIMAL &&
1868       result_status_ != MPSolver::FEASIBLE) {
1869     LOG(DFATAL) << "No solution exists. MPSolverInterface::result_status_ = "
1870                 << result_status_;
1871     return false;
1872   }
1873   return true;
1874 }
1875 
objective_value() const1876 double MPSolverInterface::objective_value() const {
1877   if (!CheckSolutionIsSynchronizedAndExists()) return 0;
1878   return objective_value_;
1879 }
1880 
best_objective_bound() const1881 double MPSolverInterface::best_objective_bound() const {
1882   const double trivial_worst_bound =
1883       maximize_ ? -std::numeric_limits<double>::infinity()
1884                 : std::numeric_limits<double>::infinity();
1885   if (!IsMIP()) {
1886     VLOG(1) << "Best objective bound only available for discrete problems.";
1887     return trivial_worst_bound;
1888   }
1889   if (!CheckSolutionIsSynchronized()) {
1890     return trivial_worst_bound;
1891   }
1892   // Special case for empty model.
1893   if (solver_->variables_.empty() && solver_->constraints_.empty()) {
1894     return solver_->Objective().offset();
1895   }
1896   return best_objective_bound_;
1897 }
1898 
InvalidateSolutionSynchronization()1899 void MPSolverInterface::InvalidateSolutionSynchronization() {
1900   if (sync_status_ == SOLUTION_SYNCHRONIZED) {
1901     sync_status_ = MODEL_SYNCHRONIZED;
1902   }
1903 }
1904 
ComputeExactConditionNumber() const1905 double MPSolverInterface::ComputeExactConditionNumber() const {
1906   // Override this method in interfaces that actually support it.
1907   LOG(DFATAL) << "ComputeExactConditionNumber not implemented for "
1908               << ProtoEnumToString<MPModelRequest::SolverType>(
1909                      static_cast<MPModelRequest::SolverType>(
1910                          solver_->ProblemType()));
1911   return 0.0;
1912 }
1913 
SetCommonParameters(const MPSolverParameters & param)1914 void MPSolverInterface::SetCommonParameters(const MPSolverParameters& param) {
1915   // TODO(user): Overhaul the code that sets parameters to enable changing
1916   // GLOP parameters without issuing warnings.
1917   // By default, we let GLOP keep its own default tolerance, much more accurate
1918   // than for the rest of the solvers.
1919   //
1920   if (solver_->ProblemType() != MPSolver::GLOP_LINEAR_PROGRAMMING) {
1921     SetPrimalTolerance(
1922         param.GetDoubleParam(MPSolverParameters::PRIMAL_TOLERANCE));
1923     SetDualTolerance(param.GetDoubleParam(MPSolverParameters::DUAL_TOLERANCE));
1924   }
1925   SetPresolveMode(param.GetIntegerParam(MPSolverParameters::PRESOLVE));
1926   // TODO(user): In the future, we could distinguish between the
1927   // algorithm to solve the root LP and the algorithm to solve node
1928   // LPs. Not sure if underlying solvers support it.
1929   int value = param.GetIntegerParam(MPSolverParameters::LP_ALGORITHM);
1930   if (value != MPSolverParameters::kDefaultIntegerParamValue) {
1931     SetLpAlgorithm(value);
1932   }
1933 }
1934 
SetMIPParameters(const MPSolverParameters & param)1935 void MPSolverInterface::SetMIPParameters(const MPSolverParameters& param) {
1936   if (solver_->ProblemType() != MPSolver::GLOP_LINEAR_PROGRAMMING) {
1937     SetRelativeMipGap(
1938         param.GetDoubleParam(MPSolverParameters::RELATIVE_MIP_GAP));
1939   }
1940 }
1941 
SetUnsupportedDoubleParam(MPSolverParameters::DoubleParam param)1942 void MPSolverInterface::SetUnsupportedDoubleParam(
1943     MPSolverParameters::DoubleParam param) {
1944   LOG(WARNING) << "Trying to set an unsupported parameter: " << param << ".";
1945 }
SetUnsupportedIntegerParam(MPSolverParameters::IntegerParam param)1946 void MPSolverInterface::SetUnsupportedIntegerParam(
1947     MPSolverParameters::IntegerParam param) {
1948   LOG(WARNING) << "Trying to set an unsupported parameter: " << param << ".";
1949 }
SetDoubleParamToUnsupportedValue(MPSolverParameters::DoubleParam param,double value)1950 void MPSolverInterface::SetDoubleParamToUnsupportedValue(
1951     MPSolverParameters::DoubleParam param, double value) {
1952   LOG(WARNING) << "Trying to set a supported parameter: " << param
1953                << " to an unsupported value: " << value;
1954 }
SetIntegerParamToUnsupportedValue(MPSolverParameters::IntegerParam param,int value)1955 void MPSolverInterface::SetIntegerParamToUnsupportedValue(
1956     MPSolverParameters::IntegerParam param, int value) {
1957   LOG(WARNING) << "Trying to set a supported parameter: " << param
1958                << " to an unsupported value: " << value;
1959 }
1960 
SetNumThreads(int num_threads)1961 absl::Status MPSolverInterface::SetNumThreads(int num_threads) {
1962   return absl::UnimplementedError(
1963       absl::StrFormat("SetNumThreads() not supported by %s.", SolverVersion()));
1964 }
1965 
SetSolverSpecificParametersAsString(const std::string & parameters)1966 bool MPSolverInterface::SetSolverSpecificParametersAsString(
1967     const std::string& parameters) {
1968   if (parameters.empty()) {
1969     return true;
1970   }
1971 
1972   LOG(WARNING) << "SetSolverSpecificParametersAsString() not supported by "
1973                << SolverVersion();
1974   return false;
1975 }
1976 
1977 // ---------- MPSolverParameters ----------
1978 
1979 const double MPSolverParameters::kDefaultRelativeMipGap = 1e-4;
1980 // For the primal and dual tolerances, choose the same default as CLP and GLPK.
1981 const double MPSolverParameters::kDefaultPrimalTolerance =
1982     operations_research::kDefaultPrimalTolerance;
1983 const double MPSolverParameters::kDefaultDualTolerance = 1e-7;
1984 const MPSolverParameters::PresolveValues MPSolverParameters::kDefaultPresolve =
1985     MPSolverParameters::PRESOLVE_ON;
1986 const MPSolverParameters::IncrementalityValues
1987     MPSolverParameters::kDefaultIncrementality =
1988         MPSolverParameters::INCREMENTALITY_ON;
1989 
1990 const double MPSolverParameters::kDefaultDoubleParamValue = -1.0;
1991 const int MPSolverParameters::kDefaultIntegerParamValue = -1;
1992 const double MPSolverParameters::kUnknownDoubleParamValue = -2.0;
1993 const int MPSolverParameters::kUnknownIntegerParamValue = -2;
1994 
1995 // The constructor sets all parameters to their default value.
MPSolverParameters()1996 MPSolverParameters::MPSolverParameters()
1997     : relative_mip_gap_value_(kDefaultRelativeMipGap),
1998       primal_tolerance_value_(kDefaultPrimalTolerance),
1999       dual_tolerance_value_(kDefaultDualTolerance),
2000       presolve_value_(kDefaultPresolve),
2001       scaling_value_(kDefaultIntegerParamValue),
2002       lp_algorithm_value_(kDefaultIntegerParamValue),
2003       incrementality_value_(kDefaultIncrementality),
2004       lp_algorithm_is_default_(true) {}
2005 
SetDoubleParam(MPSolverParameters::DoubleParam param,double value)2006 void MPSolverParameters::SetDoubleParam(MPSolverParameters::DoubleParam param,
2007                                         double value) {
2008   switch (param) {
2009     case RELATIVE_MIP_GAP: {
2010       relative_mip_gap_value_ = value;
2011       break;
2012     }
2013     case PRIMAL_TOLERANCE: {
2014       primal_tolerance_value_ = value;
2015       break;
2016     }
2017     case DUAL_TOLERANCE: {
2018       dual_tolerance_value_ = value;
2019       break;
2020     }
2021     default: {
2022       LOG(ERROR) << "Trying to set an unknown parameter: " << param << ".";
2023     }
2024   }
2025 }
2026 
SetIntegerParam(MPSolverParameters::IntegerParam param,int value)2027 void MPSolverParameters::SetIntegerParam(MPSolverParameters::IntegerParam param,
2028                                          int value) {
2029   switch (param) {
2030     case PRESOLVE: {
2031       if (value != PRESOLVE_OFF && value != PRESOLVE_ON) {
2032         LOG(ERROR) << "Trying to set a supported parameter: " << param
2033                    << " to an unknown value: " << value;
2034       }
2035       presolve_value_ = value;
2036       break;
2037     }
2038     case SCALING: {
2039       if (value != SCALING_OFF && value != SCALING_ON) {
2040         LOG(ERROR) << "Trying to set a supported parameter: " << param
2041                    << " to an unknown value: " << value;
2042       }
2043       scaling_value_ = value;
2044       break;
2045     }
2046     case LP_ALGORITHM: {
2047       if (value != DUAL && value != PRIMAL && value != BARRIER) {
2048         LOG(ERROR) << "Trying to set a supported parameter: " << param
2049                    << " to an unknown value: " << value;
2050       }
2051       lp_algorithm_value_ = value;
2052       lp_algorithm_is_default_ = false;
2053       break;
2054     }
2055     case INCREMENTALITY: {
2056       if (value != INCREMENTALITY_OFF && value != INCREMENTALITY_ON) {
2057         LOG(ERROR) << "Trying to set a supported parameter: " << param
2058                    << " to an unknown value: " << value;
2059       }
2060       incrementality_value_ = value;
2061       break;
2062     }
2063     default: {
2064       LOG(ERROR) << "Trying to set an unknown parameter: " << param << ".";
2065     }
2066   }
2067 }
2068 
ResetDoubleParam(MPSolverParameters::DoubleParam param)2069 void MPSolverParameters::ResetDoubleParam(
2070     MPSolverParameters::DoubleParam param) {
2071   switch (param) {
2072     case RELATIVE_MIP_GAP: {
2073       relative_mip_gap_value_ = kDefaultRelativeMipGap;
2074       break;
2075     }
2076     case PRIMAL_TOLERANCE: {
2077       primal_tolerance_value_ = kDefaultPrimalTolerance;
2078       break;
2079     }
2080     case DUAL_TOLERANCE: {
2081       dual_tolerance_value_ = kDefaultDualTolerance;
2082       break;
2083     }
2084     default: {
2085       LOG(ERROR) << "Trying to reset an unknown parameter: " << param << ".";
2086     }
2087   }
2088 }
2089 
ResetIntegerParam(MPSolverParameters::IntegerParam param)2090 void MPSolverParameters::ResetIntegerParam(
2091     MPSolverParameters::IntegerParam param) {
2092   switch (param) {
2093     case PRESOLVE: {
2094       presolve_value_ = kDefaultPresolve;
2095       break;
2096     }
2097     case SCALING: {
2098       scaling_value_ = kDefaultIntegerParamValue;
2099       break;
2100     }
2101     case LP_ALGORITHM: {
2102       lp_algorithm_is_default_ = true;
2103       break;
2104     }
2105     case INCREMENTALITY: {
2106       incrementality_value_ = kDefaultIncrementality;
2107       break;
2108     }
2109     default: {
2110       LOG(ERROR) << "Trying to reset an unknown parameter: " << param << ".";
2111     }
2112   }
2113 }
2114 
Reset()2115 void MPSolverParameters::Reset() {
2116   ResetDoubleParam(RELATIVE_MIP_GAP);
2117   ResetDoubleParam(PRIMAL_TOLERANCE);
2118   ResetDoubleParam(DUAL_TOLERANCE);
2119   ResetIntegerParam(PRESOLVE);
2120   ResetIntegerParam(SCALING);
2121   ResetIntegerParam(LP_ALGORITHM);
2122   ResetIntegerParam(INCREMENTALITY);
2123 }
2124 
GetDoubleParam(MPSolverParameters::DoubleParam param) const2125 double MPSolverParameters::GetDoubleParam(
2126     MPSolverParameters::DoubleParam param) const {
2127   switch (param) {
2128     case RELATIVE_MIP_GAP: {
2129       return relative_mip_gap_value_;
2130     }
2131     case PRIMAL_TOLERANCE: {
2132       return primal_tolerance_value_;
2133     }
2134     case DUAL_TOLERANCE: {
2135       return dual_tolerance_value_;
2136     }
2137     default: {
2138       LOG(ERROR) << "Trying to get an unknown parameter: " << param << ".";
2139       return kUnknownDoubleParamValue;
2140     }
2141   }
2142 }
2143 
GetIntegerParam(MPSolverParameters::IntegerParam param) const2144 int MPSolverParameters::GetIntegerParam(
2145     MPSolverParameters::IntegerParam param) const {
2146   switch (param) {
2147     case PRESOLVE: {
2148       return presolve_value_;
2149     }
2150     case LP_ALGORITHM: {
2151       if (lp_algorithm_is_default_) return kDefaultIntegerParamValue;
2152       return lp_algorithm_value_;
2153     }
2154     case INCREMENTALITY: {
2155       return incrementality_value_;
2156     }
2157     case SCALING: {
2158       return scaling_value_;
2159     }
2160     default: {
2161       LOG(ERROR) << "Trying to get an unknown parameter: " << param << ".";
2162       return kUnknownIntegerParamValue;
2163     }
2164   }
2165 }
2166 
2167 }  // namespace operations_research
2168