1 // Copyright 2010-2021 Google LLC
2 // Licensed under the Apache License, Version 2.0 (the "License");
3 // you may not use this file except in compliance with the License.
4 // You may obtain a copy of the License at
5 //
6 //     http://www.apache.org/licenses/LICENSE-2.0
7 //
8 // Unless required by applicable law or agreed to in writing, software
9 // distributed under the License is distributed on an "AS IS" BASIS,
10 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11 // See the License for the specific language governing permissions and
12 // limitations under the License.
13 
14 #include "ortools/math_opt/solvers/glop_solver.h"
15 
16 #include <algorithm>
17 #include <cstdint>
18 #include <memory>
19 #include <string>
20 #include <utility>
21 #include <vector>
22 
23 #include "ortools/base/integral_types.h"
24 #include "ortools/base/logging.h"
25 #include "absl/container/flat_hash_map.h"
26 #include "absl/memory/memory.h"
27 #include "absl/status/status.h"
28 #include "absl/status/statusor.h"
29 #include "absl/strings/str_cat.h"
30 #include "absl/strings/str_join.h"
31 #include "absl/strings/string_view.h"
32 #include "absl/time/clock.h"
33 #include "absl/time/time.h"
34 #include "absl/types/span.h"
35 #include "ortools/base/int_type.h"
36 #include "ortools/base/map_util.h"
37 #include "ortools/base/strong_vector.h"
38 #include "ortools/glop/lp_solver.h"
39 #include "ortools/glop/parameters.pb.h"
40 #include "ortools/lp_data/lp_data.h"
41 #include "ortools/lp_data/lp_types.h"
42 #include "ortools/math_opt/callback.pb.h"
43 #include "ortools/math_opt/core/math_opt_proto_utils.h"
44 #include "ortools/math_opt/core/solver_interface.h"
45 #include "ortools/math_opt/core/sparse_vector_view.h"
46 #include "ortools/math_opt/model.pb.h"
47 #include "ortools/math_opt/model_parameters.pb.h"
48 #include "ortools/math_opt/model_update.pb.h"
49 #include "ortools/math_opt/parameters.pb.h"
50 #include "ortools/math_opt/result.pb.h"
51 #include "ortools/math_opt/solution.pb.h"
52 #include "ortools/math_opt/sparse_containers.pb.h"
53 #include "ortools/port/proto_utils.h"
54 #include "absl/status/status.h"
55 #include "ortools/base/protoutil.h"
56 
57 namespace operations_research {
58 namespace math_opt {
59 
60 namespace {
61 
SafeName(const VariablesProto & variables,int index)62 absl::string_view SafeName(const VariablesProto& variables, int index) {
63   if (variables.names().empty()) {
64     return {};
65   }
66   return variables.names(index);
67 }
68 
SafeName(const LinearConstraintsProto & linear_constraints,int index)69 absl::string_view SafeName(const LinearConstraintsProto& linear_constraints,
70                            int index) {
71   if (linear_constraints.names().empty()) {
72     return {};
73   }
74   return linear_constraints.names(index);
75 }
76 
GlopVarTypeFromIsInteger(const bool is_integer)77 glop::LinearProgram::VariableType GlopVarTypeFromIsInteger(
78     const bool is_integer) {
79   return is_integer ? glop::LinearProgram::VariableType::INTEGER
80                     : glop::LinearProgram::VariableType::CONTINUOUS;
81 }
82 
83 }  // namespace
84 
GlopSolver()85 GlopSolver::GlopSolver() : linear_program_(), lp_solver_() {}
86 
AddVariables(const VariablesProto & variables)87 void GlopSolver::AddVariables(const VariablesProto& variables) {
88   for (int i = 0; i < NumVariables(variables); ++i) {
89     const glop::ColIndex col_index = linear_program_.CreateNewVariable();
90     linear_program_.SetVariableBounds(col_index, variables.lower_bounds(i),
91                                       variables.upper_bounds(i));
92     linear_program_.SetVariableName(col_index, SafeName(variables, i));
93     linear_program_.SetVariableType(
94         col_index, GlopVarTypeFromIsInteger(variables.integers(i)));
95     gtl::InsertOrDie(&variables_, variables.ids(i), col_index);
96   }
97 }
98 
99 // Note that this relies on the fact that when variable/constraint
100 // are deleted, Glop re-index everything by compacting the
101 // index domain in a stable way.
102 template <typename IndexType>
UpdateIdIndexMap(glop::StrictITIVector<IndexType,bool> indices_to_delete,IndexType num_indices,absl::flat_hash_map<int64_t,IndexType> & id_index_map)103 void UpdateIdIndexMap(glop::StrictITIVector<IndexType, bool> indices_to_delete,
104 
105                       IndexType num_indices,
106                       absl::flat_hash_map<int64_t, IndexType>& id_index_map) {
107   absl::StrongVector<IndexType, IndexType> new_indices(
108       num_indices.value(), IndexType(0));
109   IndexType new_index(0);
110   for (IndexType index(0); index < num_indices; ++index) {
111     if (indices_to_delete[index]) {
112       // Mark deleted index
113       new_indices[index] = IndexType(-1);
114     } else {
115       new_indices[index] = new_index;
116       ++new_index;
117     }
118   }
119   for (auto it = id_index_map.begin(); it != id_index_map.end();) {
120     IndexType index = it->second;
121     if (indices_to_delete[index]) {
122       // This safely deletes the entry and moves the iterator one step ahead.
123       id_index_map.erase(it++);
124     } else {
125       it->second = new_indices[index];
126       ++it;
127     }
128   }
129 }
130 
DeleteVariables(absl::Span<const int64_t> ids_to_delete)131 void GlopSolver::DeleteVariables(absl::Span<const int64_t> ids_to_delete) {
132   const glop::ColIndex num_cols = linear_program_.num_variables();
133   glop::StrictITIVector<glop::ColIndex, bool> columns_to_delete(num_cols,
134                                                                 false);
135   for (const int64_t deleted_variable_id : ids_to_delete) {
136     columns_to_delete[variables_.at(deleted_variable_id)] = true;
137   }
138   linear_program_.DeleteColumns(columns_to_delete);
139   UpdateIdIndexMap<glop::ColIndex>(columns_to_delete, num_cols, variables_);
140 }
141 
DeleteLinearConstraints(absl::Span<const int64_t> ids_to_delete)142 void GlopSolver::DeleteLinearConstraints(
143     absl::Span<const int64_t> ids_to_delete) {
144   const glop::RowIndex num_rows = linear_program_.num_constraints();
145   glop::DenseBooleanColumn rows_to_delete(num_rows, false);
146   for (const int64_t deleted_constraint_id : ids_to_delete) {
147     rows_to_delete[linear_constraints_.at(deleted_constraint_id)] = true;
148   }
149   linear_program_.DeleteRows(rows_to_delete);
150   UpdateIdIndexMap<glop::RowIndex>(rows_to_delete, num_rows,
151                                    linear_constraints_);
152 }
153 
AddLinearConstraints(const LinearConstraintsProto & linear_constraints)154 void GlopSolver::AddLinearConstraints(
155     const LinearConstraintsProto& linear_constraints) {
156   for (int i = 0; i < NumConstraints(linear_constraints); ++i) {
157     const glop::RowIndex row_index = linear_program_.CreateNewConstraint();
158     linear_program_.SetConstraintBounds(row_index,
159                                         linear_constraints.lower_bounds(i),
160                                         linear_constraints.upper_bounds(i));
161     linear_program_.SetConstraintName(row_index,
162                                       SafeName(linear_constraints, i));
163     gtl::InsertOrDie(&linear_constraints_, linear_constraints.ids(i),
164                      row_index);
165   }
166 }
167 
SetOrUpdateObjectiveCoefficients(const SparseDoubleVectorProto & linear_objective_coefficients)168 void GlopSolver::SetOrUpdateObjectiveCoefficients(
169     const SparseDoubleVectorProto& linear_objective_coefficients) {
170   for (int i = 0; i < linear_objective_coefficients.ids_size(); ++i) {
171     const glop::ColIndex col_index =
172         variables_.at(linear_objective_coefficients.ids(i));
173     linear_program_.SetObjectiveCoefficient(
174         col_index, linear_objective_coefficients.values(i));
175   }
176 }
177 
SetOrUpdateConstraintMatrix(const SparseDoubleMatrixProto & linear_constraint_matrix)178 void GlopSolver::SetOrUpdateConstraintMatrix(
179     const SparseDoubleMatrixProto& linear_constraint_matrix) {
180   for (int j = 0; j < NumMatrixNonzeros(linear_constraint_matrix); ++j) {
181     const glop::ColIndex col_index =
182         variables_.at(linear_constraint_matrix.column_ids(j));
183     const glop::RowIndex row_index =
184         linear_constraints_.at(linear_constraint_matrix.row_ids(j));
185     const double coefficient = linear_constraint_matrix.coefficients(j);
186     linear_program_.SetCoefficient(row_index, col_index, coefficient);
187   }
188 }
189 
UpdateVariableBounds(const VariableUpdatesProto & variable_updates)190 void GlopSolver::UpdateVariableBounds(
191     const VariableUpdatesProto& variable_updates) {
192   for (const auto [id, lb] : MakeView(variable_updates.lower_bounds())) {
193     const auto col_index = variables_.at(id);
194     linear_program_.SetVariableBounds(
195         col_index, lb, linear_program_.variable_upper_bounds()[col_index]);
196   }
197   for (const auto [id, ub] : MakeView(variable_updates.upper_bounds())) {
198     const auto col_index = variables_.at(id);
199     linear_program_.SetVariableBounds(
200         col_index, linear_program_.variable_lower_bounds()[col_index], ub);
201   }
202 }
203 
UpdateLinearConstraintBounds(const LinearConstraintUpdatesProto & linear_constraint_updates)204 void GlopSolver::UpdateLinearConstraintBounds(
205     const LinearConstraintUpdatesProto& linear_constraint_updates) {
206   for (const auto [id, lb] :
207        MakeView(linear_constraint_updates.lower_bounds())) {
208     const auto row_index = linear_constraints_.at(id);
209     linear_program_.SetConstraintBounds(
210         row_index, lb, linear_program_.constraint_upper_bounds()[row_index]);
211   }
212   for (const auto [id, ub] :
213        MakeView(linear_constraint_updates.upper_bounds())) {
214     const auto row_index = linear_constraints_.at(id);
215     linear_program_.SetConstraintBounds(
216         row_index, linear_program_.constraint_lower_bounds()[row_index], ub);
217   }
218 }
219 
220 std::pair<glop::GlopParameters, std::vector<std::string>>
MergeCommonParameters(const CommonSolveParametersProto & common_solver_parameters,const glop::GlopParameters & glop_parameters)221 GlopSolver::MergeCommonParameters(
222     const CommonSolveParametersProto& common_solver_parameters,
223     const glop::GlopParameters& glop_parameters) {
224   glop::GlopParameters result = glop_parameters;
225   std::vector<std::string> warnings;
226   if (!result.has_max_time_in_seconds() &&
227       common_solver_parameters.has_time_limit()) {
228     const absl::Duration time_limit =
229         util_time::DecodeGoogleApiProto(common_solver_parameters.time_limit())
230             .value();
231     result.set_max_time_in_seconds(absl::ToDoubleSeconds(time_limit));
232   }
233   if (!result.has_log_search_progress()) {
234     result.set_log_search_progress(common_solver_parameters.enable_output());
235   }
236   if (!result.has_num_omp_threads() && common_solver_parameters.has_threads()) {
237     result.set_num_omp_threads(common_solver_parameters.threads());
238   }
239   if (!result.has_random_seed() && common_solver_parameters.has_random_seed()) {
240     const int random_seed = std::max(0, common_solver_parameters.random_seed());
241     result.set_random_seed(random_seed);
242   }
243   if (!result.has_use_dual_simplex() &&
244       common_solver_parameters.lp_algorithm() != LP_ALGORITHM_UNSPECIFIED) {
245     switch (common_solver_parameters.lp_algorithm()) {
246       case LP_ALGORITHM_PRIMAL_SIMPLEX:
247         result.set_use_dual_simplex(false);
248         break;
249       case LP_ALGORITHM_DUAL_SIMPLEX:
250         result.set_use_dual_simplex(true);
251         break;
252       case LP_ALGORITHM_BARRIER:
253         warnings.push_back(
254             "GLOP does not support 'LP_ALGORITHM_BARRIER' value for "
255             "'lp_algorithm' parameter.");
256         break;
257       default:
258         LOG(FATAL) << "LPAlgorithm: "
259                    << ProtoEnumToString(common_solver_parameters.lp_algorithm())
260                    << " unknown, error setting GLOP parameters";
261     }
262   }
263   if (!result.has_use_scaling() && !result.has_scaling_method() &&
264       common_solver_parameters.scaling() != EMPHASIS_UNSPECIFIED) {
265     switch (common_solver_parameters.scaling()) {
266       case EMPHASIS_OFF:
267         result.set_use_scaling(false);
268         break;
269       case EMPHASIS_LOW:
270       case EMPHASIS_MEDIUM:
271         result.set_use_scaling(true);
272         result.set_scaling_method(glop::GlopParameters::EQUILIBRATION);
273         break;
274       case EMPHASIS_HIGH:
275       case EMPHASIS_VERY_HIGH:
276         result.set_use_scaling(true);
277         result.set_scaling_method(glop::GlopParameters::LINEAR_PROGRAM);
278         break;
279       default:
280         LOG(FATAL) << "Scaling emphasis: "
281                    << ProtoEnumToString(common_solver_parameters.scaling())
282                    << " unknown, error setting GLOP parameters";
283     }
284   }
285   if (!result.has_use_preprocessing() &&
286       common_solver_parameters.presolve() != EMPHASIS_UNSPECIFIED) {
287     switch (common_solver_parameters.presolve()) {
288       case EMPHASIS_OFF:
289         result.set_use_preprocessing(false);
290         break;
291       case EMPHASIS_LOW:
292       case EMPHASIS_MEDIUM:
293       case EMPHASIS_HIGH:
294       case EMPHASIS_VERY_HIGH:
295         result.set_use_preprocessing(true);
296         break;
297       default:
298         LOG(FATAL) << "Presolve emphasis: "
299                    << ProtoEnumToString(common_solver_parameters.presolve())
300                    << " unknown, error setting GLOP parameters";
301     }
302   }
303   if (common_solver_parameters.cuts() != EMPHASIS_UNSPECIFIED) {
304     warnings.push_back(absl::StrCat(
305         "GLOP does not support 'cuts' parameters, but cuts was set to: ",
306         ProtoEnumToString(common_solver_parameters.cuts())));
307   }
308   if (common_solver_parameters.heuristics() != EMPHASIS_UNSPECIFIED) {
309     warnings.push_back(
310         absl::StrCat("GLOP does not support 'heuristics' parameter, but "
311                      "heuristics was set to: ",
312                      ProtoEnumToString(common_solver_parameters.heuristics())));
313   }
314   return std::make_pair(std::move(result), std::move(warnings));
315 }
316 
CanUpdate(const ModelUpdateProto & model_update)317 bool GlopSolver::CanUpdate(const ModelUpdateProto& model_update) {
318   return true;
319 }
320 
321 template <typename IndexType>
FillSparseDoubleVector(const std::vector<int64_t> & ids_in_order,const absl::flat_hash_map<int64_t,IndexType> & id_map,const glop::StrictITIVector<IndexType,glop::Fractional> & values,const SparseVectorFilterProto & filter)322 SparseDoubleVectorProto FillSparseDoubleVector(
323     const std::vector<int64_t>& ids_in_order,
324     const absl::flat_hash_map<int64_t, IndexType>& id_map,
325     const glop::StrictITIVector<IndexType, glop::Fractional>& values,
326     const SparseVectorFilterProto& filter) {
327   SparseVectorFilterPredicate predicate(filter);
328   SparseDoubleVectorProto result;
329   for (const int64_t variable_id : ids_in_order) {
330     const double value = values[id_map.at(variable_id)];
331     if (predicate.AcceptsAndUpdate(variable_id, value)) {
332       result.add_ids(variable_id);
333       result.add_values(value);
334     }
335   }
336   return result;
337 }
338 
339 template <typename T>
GetSortedIs(const absl::flat_hash_map<int64_t,T> & id_map)340 std::vector<int64_t> GetSortedIs(
341     const absl::flat_hash_map<int64_t, T>& id_map) {
342   std::vector<int64_t> sorted;
343   sorted.reserve(id_map.size());
344   for (const auto& entry : id_map) {
345     sorted.emplace_back(entry.first);
346   }
347   std::sort(sorted.begin(), sorted.end());
348   return sorted;
349 }
350 
FillSolveResult(const glop::ProblemStatus status,const ModelSolveParametersProto & model_parameters,SolveResultProto & solve_result)351 void GlopSolver::FillSolveResult(
352     const glop::ProblemStatus status,
353     const ModelSolveParametersProto& model_parameters,
354     SolveResultProto& solve_result) {
355   solve_result.mutable_solve_stats()->set_simplex_iterations(
356       lp_solver_.GetNumberOfSimplexIterations());
357   // TODO(b/168374742): this needs to be properly filled in. In particular, we
358   // can give better primal and dual bounds when the status is not OPTIMAL.
359   const bool is_maximize = linear_program_.IsMaximizationProblem();
360   constexpr double kInf = std::numeric_limits<double>::infinity();
361   solve_result.mutable_solve_stats()->set_best_primal_bound(is_maximize ? -kInf
362                                                                         : kInf);
363   solve_result.mutable_solve_stats()->set_best_dual_bound(is_maximize ? kInf
364                                                                       : -kInf);
365   if (status == glop::ProblemStatus::OPTIMAL) {
366     solve_result.set_termination_reason(SolveResultProto::OPTIMAL);
367     solve_result.mutable_solve_stats()->set_best_primal_bound(
368         lp_solver_.GetObjectiveValue());
369     solve_result.mutable_solve_stats()->set_best_dual_bound(
370         lp_solver_.GetObjectiveValue());
371 
372     auto sorted_variables = GetSortedIs(variables_);
373     auto sorted_constraints = GetSortedIs(linear_constraints_);
374 
375     PrimalSolutionProto* const primal_solution =
376         solve_result.add_primal_solutions();
377     primal_solution->set_objective_value(lp_solver_.GetObjectiveValue());
378     *primal_solution->mutable_variable_values() = FillSparseDoubleVector(
379         sorted_variables, variables_, lp_solver_.variable_values(),
380         model_parameters.primal_variables_filter());
381 
382     DualSolutionProto* const dual_solution = solve_result.add_dual_solutions();
383     dual_solution->set_objective_value(lp_solver_.GetObjectiveValue());
384     *dual_solution->mutable_dual_values() = FillSparseDoubleVector(
385         sorted_constraints, linear_constraints_, lp_solver_.dual_values(),
386         model_parameters.dual_linear_constraints_filter());
387     *dual_solution->mutable_reduced_costs() = FillSparseDoubleVector(
388         sorted_variables, variables_, lp_solver_.reduced_costs(),
389         model_parameters.dual_variables_filter());
390     // TODO(user): consider pulling these out to a separate method once we
391     // support all statuses
392   } else if (status == glop::ProblemStatus::PRIMAL_INFEASIBLE ||
393              status == glop::ProblemStatus::DUAL_UNBOUNDED) {
394     solve_result.set_termination_reason(SolveResultProto::INFEASIBLE);
395   } else if (status == glop::ProblemStatus::PRIMAL_UNBOUNDED) {
396     solve_result.set_termination_reason(SolveResultProto::UNBOUNDED);
397   } else if (status == glop::ProblemStatus::DUAL_INFEASIBLE ||
398              status == glop::ProblemStatus::INFEASIBLE_OR_UNBOUNDED) {
399     solve_result.set_termination_reason(SolveResultProto::DUAL_INFEASIBLE);
400   } else {
401     LOG(DFATAL) << "Termination not implemented.";
402     solve_result.set_termination_reason(
403         SolveResultProto::TERMINATION_REASON_UNSPECIFIED);
404     solve_result.set_termination_detail(absl::StrCat("Glop status: ", status));
405   }
406 }
407 
Solve(const SolveParametersProto & parameters,const ModelSolveParametersProto & model_parameters,const CallbackRegistrationProto & callback_registration,const Callback cb)408 absl::StatusOr<SolveResultProto> GlopSolver::Solve(
409     const SolveParametersProto& parameters,
410     const ModelSolveParametersProto& model_parameters,
411     const CallbackRegistrationProto& callback_registration, const Callback cb) {
412   const absl::Time start = absl::Now();
413   SolveResultProto result;
414   {
415     auto [glop_parameters, warnings] = MergeCommonParameters(
416         parameters.common_parameters(), parameters.glop_parameters());
417     if (!warnings.empty()) {
418       if (parameters.common_parameters().strictness().bad_parameter()) {
419         return absl::InvalidArgumentError(absl::StrJoin(warnings, "; "));
420       } else {
421         for (std::string& warning : warnings) {
422           result.add_warnings(std::move(warning));
423         }
424       }
425     }
426     lp_solver_.SetParameters(glop_parameters);
427   }
428 
429   const glop::ProblemStatus status = lp_solver_.Solve(linear_program_);
430   FillSolveResult(status, model_parameters, result);
431   CHECK_OK(util_time::EncodeGoogleApiProto(
432       absl::Now() - start, result.mutable_solve_stats()->mutable_solve_time()));
433   return result;
434 }
435 
New(const ModelProto & model,const SolverInitializerProto & initializer)436 absl::StatusOr<std::unique_ptr<SolverInterface>> GlopSolver::New(
437     const ModelProto& model, const SolverInitializerProto& initializer) {
438   auto solver = absl::WrapUnique(new GlopSolver);
439   solver->linear_program_.SetName(model.name());
440   solver->linear_program_.SetMaximizationProblem(model.objective().maximize());
441   solver->linear_program_.SetObjectiveOffset(model.objective().offset());
442 
443   solver->AddVariables(model.variables());
444   solver->SetOrUpdateObjectiveCoefficients(
445       model.objective().linear_coefficients());
446 
447   solver->AddLinearConstraints(model.linear_constraints());
448   solver->SetOrUpdateConstraintMatrix(model.linear_constraint_matrix());
449   solver->linear_program_.CleanUp();
450   return solver;
451 }
452 
Update(const ModelUpdateProto & model_update)453 absl::Status GlopSolver::Update(const ModelUpdateProto& model_update) {
454   if (model_update.objective_updates().has_direction_update()) {
455     linear_program_.SetMaximizationProblem(
456         model_update.objective_updates().direction_update());
457   }
458   if (model_update.objective_updates().has_offset_update()) {
459     linear_program_.SetObjectiveOffset(
460         model_update.objective_updates().offset_update());
461   }
462 
463   DeleteVariables(model_update.deleted_variable_ids());
464   AddVariables(model_update.new_variables());
465 
466   SetOrUpdateObjectiveCoefficients(
467       model_update.objective_updates().linear_coefficients());
468   UpdateVariableBounds(model_update.variable_updates());
469 
470   DeleteLinearConstraints(model_update.deleted_linear_constraint_ids());
471   AddLinearConstraints(model_update.new_linear_constraints());
472   UpdateLinearConstraintBounds(model_update.linear_constraint_updates());
473 
474   SetOrUpdateConstraintMatrix(model_update.linear_constraint_matrix_updates());
475 
476   linear_program_.CleanUp();
477 
478   return absl::OkStatus();
479 }
480 
481 MATH_OPT_REGISTER_SOLVER(SOLVER_TYPE_GLOP, GlopSolver::New)
482 
483 }  // namespace math_opt
484 }  // namespace operations_research
485