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