From 558cbf9cd447ac89c8858424c55ff01bb689d8bb Mon Sep 17 00:00:00 2001 From: Corentin Le Molgat Date: Thu, 31 Mar 2022 15:30:39 +0200 Subject: [PATCH] math_opt: sync --- ortools/math_opt/core/math_opt_proto_utils.cc | 1 + ortools/math_opt/core/model_storage.cc | 9 +- ortools/math_opt/core/model_storage.h | 5 +- ortools/math_opt/core/solver.cc | 26 ++ ortools/math_opt/core/solver.h | 5 + ortools/math_opt/core/solver_interface.h | 5 +- ortools/math_opt/cpp/BUILD.bazel | 10 + ortools/math_opt/cpp/linear_constraint.h | 8 +- ortools/math_opt/cpp/matchers.cc | 45 ++++ ortools/math_opt/cpp/matchers.h | 52 +++- ortools/math_opt/cpp/model.cc | 5 +- ortools/math_opt/cpp/model.h | 5 +- ortools/math_opt/cpp/solve.cc | 10 +- ortools/math_opt/cpp/statistics.cc | 117 +++++++++ ortools/math_opt/cpp/statistics.h | 66 +++++ .../math_opt/cpp/variable_and_expressions.h | 245 ++++++++++++++++-- ortools/math_opt/io/BUILD.bazel | 10 + ortools/math_opt/io/names_removal.cc | 32 +++ ortools/math_opt/io/names_removal.h | 36 +++ ortools/math_opt/parameters.proto | 2 + ortools/math_opt/samples/basic_example.cc | 4 +- ortools/math_opt/samples/cocktail_hour.cc | 2 +- .../math_opt/samples/facility_lp_benders.cc | 4 +- .../math_opt/samples/integer_programming.cc | 2 +- ortools/math_opt/samples/tsp.cc | 3 + ortools/math_opt/solvers/BUILD.bazel | 1 + ortools/math_opt/solvers/glpk/BUILD.bazel | 7 + .../solvers/glpk/glpk_sparse_vector.cc | 73 ++++++ .../solvers/glpk/glpk_sparse_vector.h | 201 ++++++++++++++ ortools/math_opt/solvers/glpk/rays.cc | 4 +- ortools/math_opt/solvers/glpk_solver.cc | 143 ++-------- ortools/math_opt/solvers/gurobi/g_gurobi.cc | 2 +- ortools/math_opt/solvers/gurobi/g_gurobi.h | 3 +- ortools/math_opt/solvers/gurobi_solver.cc | 48 ++-- ortools/math_opt/solvers/gurobi_solver.h | 1 + ortools/math_opt/solvers/pdlp_solver.h | 3 - ortools/math_opt/tools/BUILD.bazel | 2 + ortools/math_opt/tools/mathopt_solve_main.cc | 100 ++++++- .../math_opt/validators/result_validator.cc | 15 +- 39 files changed, 1105 insertions(+), 207 deletions(-) create mode 100644 ortools/math_opt/cpp/statistics.cc create mode 100644 ortools/math_opt/cpp/statistics.h create mode 100644 ortools/math_opt/io/names_removal.cc create mode 100644 ortools/math_opt/io/names_removal.h create mode 100644 ortools/math_opt/solvers/glpk/glpk_sparse_vector.cc create mode 100644 ortools/math_opt/solvers/glpk/glpk_sparse_vector.h diff --git a/ortools/math_opt/core/math_opt_proto_utils.cc b/ortools/math_opt/core/math_opt_proto_utils.cc index 917b766191..56c45c1830 100644 --- a/ortools/math_opt/core/math_opt_proto_utils.cc +++ b/ortools/math_opt/core/math_opt_proto_utils.cc @@ -17,6 +17,7 @@ #include #include +#include #include "absl/container/flat_hash_set.h" #include "ortools/base/integral_types.h" diff --git a/ortools/math_opt/core/model_storage.cc b/ortools/math_opt/core/model_storage.cc index 1b2d534ae8..5409fb274a 100644 --- a/ortools/math_opt/core/model_storage.cc +++ b/ortools/math_opt/core/model_storage.cc @@ -204,9 +204,14 @@ void ModelStorage::UpdateLinearConstraintCoefficients( } } -std::unique_ptr ModelStorage::Clone() const { +std::unique_ptr ModelStorage::Clone( + const std::optional new_name) const { + ModelProto model_proto = ExportModel(); + if (new_name.has_value()) { + model_proto.set_name(std::string(*new_name)); + } absl::StatusOr> clone = - ModelStorage::FromModelProto(ExportModel()); + ModelStorage::FromModelProto(model_proto); // Unless there is a very serious bug, a model exported by ExportModel() // should always be valid. CHECK_OK(clone.status()); diff --git a/ortools/math_opt/core/model_storage.h b/ortools/math_opt/core/model_storage.h index 271f6c62d9..f7b2ffca16 100644 --- a/ortools/math_opt/core/model_storage.h +++ b/ortools/math_opt/core/model_storage.h @@ -170,13 +170,14 @@ class ModelStorage { ModelStorage(const ModelStorage&) = delete; ModelStorage& operator=(const ModelStorage&) = delete; - // Returns a clone of the model. + // Returns a clone of the model, optionally changing model's name. // // The variables and constraints have the same ids. The clone will also not // reused any id of variable/constraint that was deleted in the original. // // Note that the returned model does not have any update tracker. - std::unique_ptr Clone() const; + std::unique_ptr Clone( + std::optional new_name = std::nullopt) const; inline const std::string& name() const { return name_; } diff --git a/ortools/math_opt/core/solver.cc b/ortools/math_opt/core/solver.cc index 8fa56912f1..6e9bc14dda 100644 --- a/ortools/math_opt/core/solver.cc +++ b/ortools/math_opt/core/solver.cc @@ -156,6 +156,14 @@ class ConcurrentCallsGuard { absl::Mutex* mutex_; }; +// Returns the Status returned by Solve() & Update() when called after a +// previous call to one of them failed. +absl::Status PreviousFatalFailureOccurred() { + return absl::InvalidArgumentError( + "a previous call to Solve() or Update() failed, the Solver can't be used " + "anymore"); +} + } // namespace absl::StatusOr Solver::NonIncrementalSolve( @@ -192,6 +200,13 @@ absl::StatusOr> Solver::New( absl::StatusOr Solver::Solve(const SolveArgs& arguments) { ASSIGN_OR_RETURN(const auto guard, ConcurrentCallsGuard::TryAcquire(mutex_)); + if (fatal_failure_occurred_) { + return PreviousFatalFailureOccurred(); + } + + // We will reset it in code paths where no error occur. + fatal_failure_occurred_ = true; + // TODO(b/168037341): we should validate the result maths. Since the result // can be filtered, this should be included in the solver_interface // implementations. @@ -231,18 +246,29 @@ absl::StatusOr Solver::Solve(const SolveArgs& arguments) { RETURN_IF_ERROR(ToInternalError( ValidateResult(result, arguments.model_parameters, model_summary_))); + fatal_failure_occurred_ = false; return result; } absl::StatusOr Solver::Update(const ModelUpdateProto& model_update) { ASSIGN_OR_RETURN(const auto guard, ConcurrentCallsGuard::TryAcquire(mutex_)); + if (fatal_failure_occurred_) { + return PreviousFatalFailureOccurred(); + } + + // We will reset it in code paths where no error occur. + fatal_failure_occurred_ = true; + RETURN_IF_ERROR(ValidateModelUpdateAndSummary(model_update, model_summary_)); if (!underlying_solver_->CanUpdate(model_update)) { + fatal_failure_occurred_ = false; return false; } UpdateSummaryFromModelUpdate(model_update, model_summary_); RETURN_IF_ERROR(underlying_solver_->Update(model_update)); + + fatal_failure_occurred_ = false; return true; } diff --git a/ortools/math_opt/core/solver.h b/ortools/math_opt/core/solver.h index 0675fde220..b0e7612f5b 100644 --- a/ortools/math_opt/core/solver.h +++ b/ortools/math_opt/core/solver.h @@ -135,7 +135,12 @@ class Solver { absl::Mutex mutex_; const std::unique_ptr underlying_solver_; + ModelSummary model_summary_; + + // Set to true if a previous call to Solve() or Update() returned a failing + // status. + bool fatal_failure_occurred_ = false; }; namespace internal { diff --git a/ortools/math_opt/core/solver_interface.h b/ortools/math_opt/core/solver_interface.h index 91a18eedfe..2dcba61767 100644 --- a/ortools/math_opt/core/solver_interface.h +++ b/ortools/math_opt/core/solver_interface.h @@ -49,8 +49,9 @@ inline constexpr absl::string_view kMessageCallbackNotSupported = // // This interface is not meant to be used directly. The actual API is the one of // the Solver class. The Solver class validates the models before calling this -// interface. It also makes sure no concurrent calls happen on Solve(), -// CanUpdate() and Update(). +// interface. It makes sure no concurrent calls happen on Solve(), CanUpdate() +// and Update(). It makes sure no other function is called after Solve(), +// Update() or a callback have failed. // // Implementations of this interface should not have public constructors but // instead have a static `New` function with the signature of Factory function diff --git a/ortools/math_opt/cpp/BUILD.bazel b/ortools/math_opt/cpp/BUILD.bazel index c2af36300c..0c14241f4e 100644 --- a/ortools/math_opt/cpp/BUILD.bazel +++ b/ortools/math_opt/cpp/BUILD.bazel @@ -326,3 +326,13 @@ cc_library( "@com_google_absl//absl/types:span", ], ) + +cc_library( + name = "statistics", + srcs = ["statistics.cc"], + hdrs = ["statistics.h"], + deps = [ + ":model", + "//ortools/math_opt/core:model_storage", + ], +) diff --git a/ortools/math_opt/cpp/linear_constraint.h b/ortools/math_opt/cpp/linear_constraint.h index d112f1a811..342d0e1031 100644 --- a/ortools/math_opt/cpp/linear_constraint.h +++ b/ortools/math_opt/cpp/linear_constraint.h @@ -124,7 +124,13 @@ H AbslHashValue(H h, const LinearConstraint& linear_constraint) { std::ostream& operator<<(std::ostream& ostr, const LinearConstraint& linear_constraint) { - ostr << linear_constraint.name(); + // TODO(b/170992529): handle quoting of invalid characters in the name. + const std::string& name = linear_constraint.name(); + if (name.empty()) { + ostr << "__lin_con#" << linear_constraint.id() << "__"; + } else { + ostr << name; + } return ostr; } diff --git a/ortools/math_opt/cpp/matchers.cc b/ortools/math_opt/cpp/matchers.cc index 1d7f34e0fc..321a168be3 100644 --- a/ortools/math_opt/cpp/matchers.cc +++ b/ortools/math_opt/cpp/matchers.cc @@ -28,6 +28,7 @@ #include "gtest/gtest.h" #include "ortools/base/logging.h" #include "ortools/math_opt/cpp/math_opt.h" +#include "ortools/math_opt/cpp/variable_and_expressions.h" namespace operations_research { namespace math_opt { @@ -40,6 +41,7 @@ using ::testing::AnyOf; using ::testing::AnyOfArray; using ::testing::Contains; using ::testing::DoubleNear; +using ::testing::Eq; using ::testing::ExplainMatchResult; using ::testing::Field; using ::testing::IsEmpty; @@ -48,6 +50,7 @@ using ::testing::MatcherInterface; using ::testing::MatchResultListener; using ::testing::Optional; using ::testing::PrintToString; +using ::testing::Property; } // namespace //////////////////////////////////////////////////////////////////////////////// @@ -231,6 +234,48 @@ Matcher> IsNear( tolerance)); } +template +Matcher> IsNear(IdMap expected, + const double tolerance) { + return Matcher>( + new IdMapMatcher(std::move(expected), /*all_keys=*/true, tolerance)); +} + +template +Matcher> IsNearlySubsetOf(IdMap expected, + const double tolerance) { + return Matcher>( + new IdMapMatcher(std::move(expected), /*all_keys=*/false, tolerance)); +} + +//////////////////////////////////////////////////////////////////////////////// +// Matchers for LinearExpression and QuadraticExpression +//////////////////////////////////////////////////////////////////////////////// + +testing::Matcher IsIdentical(LinearExpression expected) { + CHECK(!std::isnan(expected.offset())) << "Illegal NaN-valued offset"; + return AllOf( + Property("storage", &LinearExpression::storage, Eq(expected.storage())), + Property("offset", &LinearExpression::offset, + testing::Eq(expected.offset())), + Property("terms", &LinearExpression::terms, + IsNear(expected.terms(), /*tolerance=*/0))); +} + +testing::Matcher IsIdentical( + QuadraticExpression expected) { + CHECK(!std::isnan(expected.offset())) << "Illegal NaN-valued offset"; + return AllOf( + Property("storage", &QuadraticExpression::storage, + Eq(expected.storage())), + Property("offset", &QuadraticExpression::offset, + testing::Eq(expected.offset())), + Property("linear_terms", &QuadraticExpression::linear_terms, + IsNear(expected.linear_terms(), /*tolerance=*/0)), + Property("quadratic_terms", &QuadraticExpression::quadratic_terms, + IsNear(expected.quadratic_terms(), /*tolerance=*/0))); +} + //////////////////////////////////////////////////////////////////////////////// // Matcher helpers //////////////////////////////////////////////////////////////////////////////// diff --git a/ortools/math_opt/cpp/matchers.h b/ortools/math_opt/cpp/matchers.h index 26606f32a1..40ce4b9bb4 100644 --- a/ortools/math_opt/cpp/matchers.h +++ b/ortools/math_opt/cpp/matchers.h @@ -106,29 +106,71 @@ namespace math_opt { constexpr double kMatcherDefaultTolerance = 1e-5; //////////////////////////////////////////////////////////////////////////////// -// Matchers for IdMap and IdMap +// Matchers for IdMap //////////////////////////////////////////////////////////////////////////////// -// Checks that the maps have identical keys and values within tolerance. +// Checks that the maps have identical keys and values within tolerance. This +// factory will CHECK-fail if expected contains any NaN values. testing::Matcher> IsNear( VariableMap expected, double tolerance = kMatcherDefaultTolerance); // Checks that the keys of actual are a subset of the keys of expected, and that -// for all shared keys, the values are within tolerance. +// for all shared keys, the values are within tolerance. This factory will +// CHECK-fail if expected contains any NaN values, and any NaN values in the +// expression compared against will result in the matcher failing. testing::Matcher> IsNearlySubsetOf( VariableMap expected, double tolerance = kMatcherDefaultTolerance); -// Checks that the maps have identical keys and values within tolerance. +// Checks that the maps have identical keys and values within tolerance. This +// factory will CHECK-fail if expected contains any NaN values, and any NaN +// values in the expression compared against will result in the matcher failing. testing::Matcher> IsNear( LinearConstraintMap expected, double tolerance = kMatcherDefaultTolerance); // Checks that the keys of actual are a subset of the keys of expected, and that -// for all shared keys, the values are within tolerance. +// for all shared keys, the values are within tolerance. This factory will +// CHECK-fail if expected contains any NaN values, and any NaN values in the +// expression compared against will result in the matcher failing. testing::Matcher> IsNearlySubsetOf( LinearConstraintMap expected, double tolerance = kMatcherDefaultTolerance); +// Checks that the maps have identical keys and values within tolerance. Works +// for VariableMap, LinearConstraintMap, among other realizations of IdMap. This +// factory will CHECK-fail if expected contains any NaN values, and any NaN +// values in the expression compared against will result in the matcher failing. +template +testing::Matcher> IsNear( + IdMap expected, + const double tolerance = kMatcherDefaultTolerance); + +// Checks that the keys of actual are a subset of the keys of expected, and that +// for all shared keys, the values are within tolerance. Works for VariableMap, +// LinearConstraintMap, among other realizations of IdMap. This factory will +// CHECK-fail if expected contains any NaN values, and any NaN values in the +// expression compared against will result in the matcher failing. +template +testing::Matcher> IsNearlySubsetOf( + IdMap expected, + const double tolerance = kMatcherDefaultTolerance); + +//////////////////////////////////////////////////////////////////////////////// +// Matchers for LinearExpression and QuadraticExpression +//////////////////////////////////////////////////////////////////////////////// + +// Checks that the expressions are structurally identical (i.e., internal maps +// have the same keys and storage, coefficients are exactly equal). This factory +// will CHECK-fail if expected contains any NaN values, and any NaN values in +// the expression compared against will result in the matcher failing. +testing::Matcher IsIdentical(LinearExpression expected); + +// Checks that the expressions are structurally identical (i.e., internal maps +// have the same keys and storage, coefficients are exactly equal). This factory +// will CHECK-fail if expected contains any NaN values, and any NaN values in +// the expression compared against will result in the matcher failing. +testing::Matcher IsIdentical(QuadraticExpression expected); + //////////////////////////////////////////////////////////////////////////////// // Matchers for solutions //////////////////////////////////////////////////////////////////////////////// diff --git a/ortools/math_opt/cpp/model.cc b/ortools/math_opt/cpp/model.cc index 5990e5db5d..36fb12fff7 100644 --- a/ortools/math_opt/cpp/model.cc +++ b/ortools/math_opt/cpp/model.cc @@ -44,8 +44,9 @@ Model::Model(const absl::string_view name) Model::Model(std::unique_ptr storage) : storage_(std::move(storage)) {} -std::unique_ptr Model::Clone() const { - return std::make_unique(storage_->Clone()); +std::unique_ptr Model::Clone( + const std::optional new_name) const { + return std::make_unique(storage_->Clone(new_name)); } LinearConstraint Model::AddLinearConstraint( diff --git a/ortools/math_opt/cpp/model.h b/ortools/math_opt/cpp/model.h index c0be707646..587e130b83 100644 --- a/ortools/math_opt/cpp/model.h +++ b/ortools/math_opt/cpp/model.h @@ -120,7 +120,7 @@ class Model { Model(const Model&) = delete; Model& operator=(const Model&) = delete; - // Returns a clone of this model. + // Returns a clone of this model, optionally changing the model's name. // // The variables and constraints have the same integer ids. The clone will // also not reused any id of variable/constraint that was deleted in the @@ -137,7 +137,8 @@ class Model { // * in an arbitrary order using Variables() and LinearConstraints(). // // Note that the returned model does not have any update tracker. - std::unique_ptr Clone() const; + std::unique_ptr Clone( + std::optional new_name = std::nullopt) const; inline const std::string& name() const; diff --git a/ortools/math_opt/cpp/solve.cc b/ortools/math_opt/cpp/solve.cc index 0e87f60c2c..40825f37dd 100644 --- a/ortools/math_opt/cpp/solve.cc +++ b/ortools/math_opt/cpp/solve.cc @@ -27,6 +27,7 @@ #include "ortools/math_opt/core/model_storage.h" #include "ortools/math_opt/core/solver.h" #include "ortools/math_opt/cpp/model.h" +#include "ortools/util/status_macros.h" namespace operations_research { namespace math_opt { @@ -141,7 +142,8 @@ absl::StatusOr IncrementalSolver::Update() { return UpdateResult(true, std::move(model_update)); } - ASSIGN_OR_RETURN(const bool did_update, solver_->Update(*model_update)); + OR_ASSIGN_OR_RETURN3(const bool did_update, solver_->Update(*model_update), + _ << "update failed"); RETURN_IF_ERROR(update_tracker_->Checkpoint()); if (did_update) { @@ -150,8 +152,10 @@ absl::StatusOr IncrementalSolver::Update() { ASSIGN_OR_RETURN(const ModelProto model_proto, update_tracker_->ExportModel()); - ASSIGN_OR_RETURN(solver_, Solver::New(EnumToProto(solver_type_), model_proto, - ToSolverInitArgs(init_args_))); + OR_ASSIGN_OR_RETURN3(solver_, + Solver::New(EnumToProto(solver_type_), model_proto, + ToSolverInitArgs(init_args_)), + _ << "solver re-creation failed"); return UpdateResult(false, std::move(model_update)); } diff --git a/ortools/math_opt/cpp/statistics.cc b/ortools/math_opt/cpp/statistics.cc new file mode 100644 index 0000000000..2ffca08ec5 --- /dev/null +++ b/ortools/math_opt/cpp/statistics.cc @@ -0,0 +1,117 @@ +// Copyright 2010-2021 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "ortools/math_opt/cpp/statistics.h" + +#include +#include +#include +#include +#include +#include + +#include "ortools/math_opt/core/model_storage.h" +#include "ortools/math_opt/cpp/model.h" + +namespace operations_research::math_opt { +namespace { + +// Formatter for std::optional that uses the given width for std::setw(). +struct OptionalRangeFormatter { + OptionalRangeFormatter(const std::optional& range, const int width) + : range(range), width(width) {} + + const std::optional& range; + const int width; +}; + +std::ostream& operator<<(std::ostream& out, const OptionalRangeFormatter& fmt) { + if (!fmt.range.has_value()) { + out << "no finite values"; + return out; + } + + out << '[' << std::setw(fmt.width) << fmt.range->first << ", " + << std::setw(fmt.width) << fmt.range->second << ']'; + + return out; +} + +// Updates the input optional range with abs(v) if it is finite and non-zero. +void UpdateOptionalRange(std::optional& range, const double v) { + if (std::isinf(v) || v == 0.0) { + return; + } + + const double abs_v = std::abs(v); + if (range.has_value()) { + range->first = std::min(range->first, abs_v); + range->second = std::max(range->second, abs_v); + } else { + range = std::make_pair(abs_v, abs_v); + } +} + +} // namespace + +std::ostream& operator<<(std::ostream& out, const ModelRanges& ranges) { + const auto last_precision = out.precision(2); + const auto last_flags = out.flags(); + out.setf(std::ios_base::scientific, std::ios_base::floatfield); + out.setf(std::ios_base::left, std::ios_base::adjustfield); + // Numbers are printed in scientific notation with a precision of 2. Since + // they are expected to be positive we can ignore the optional leading minus + // sign. We thus expects `d.dde[+-]dd(d)?` (the exponent is at least 2 digits + // but double can require 3 digits, with max +308 and min -308). Thus we can + // use a width of 9 to align the ranges properly. + constexpr int kWidth = 9; + + out << "Objective terms : " + << OptionalRangeFormatter(ranges.objective_terms, kWidth) + << "\nVariable bounds : " + << OptionalRangeFormatter(ranges.variable_bounds, kWidth) + << "\nLinear constraints bounds : " + << OptionalRangeFormatter(ranges.linear_constraint_bounds, kWidth) + << "\nLinear constraints coeffs : " + << OptionalRangeFormatter(ranges.linear_constraint_coefficients, kWidth); + + out.precision(last_precision); + out.flags(last_flags); + + return out; +} + +ModelRanges ComputeModelRanges(const Model& model) { + ModelRanges ranges; + const auto objective = model.ObjectiveAsQuadraticExpression(); + for (const auto& [_, coeff] : objective.linear_terms()) { + UpdateOptionalRange(ranges.objective_terms, coeff); + } + for (const auto& [_, coeff] : objective.quadratic_terms()) { + UpdateOptionalRange(ranges.objective_terms, coeff); + } + for (const Variable& v : model.Variables()) { + UpdateOptionalRange(ranges.variable_bounds, v.lower_bound()); + UpdateOptionalRange(ranges.variable_bounds, v.upper_bound()); + } + for (const LinearConstraint& c : model.LinearConstraints()) { + UpdateOptionalRange(ranges.linear_constraint_bounds, c.lower_bound()); + UpdateOptionalRange(ranges.linear_constraint_bounds, c.upper_bound()); + } + for (const auto& [_, coeff] : model.storage()->linear_constraint_matrix()) { + UpdateOptionalRange(ranges.linear_constraint_coefficients, coeff); + } + return ranges; +} + +} // namespace operations_research::math_opt diff --git a/ortools/math_opt/cpp/statistics.h b/ortools/math_opt/cpp/statistics.h new file mode 100644 index 0000000000..665f8a6cb7 --- /dev/null +++ b/ortools/math_opt/cpp/statistics.h @@ -0,0 +1,66 @@ +// Copyright 2010-2021 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#ifndef OR_TOOLS_MATH_OPT_CPP_STATISTICS_H_ +#define OR_TOOLS_MATH_OPT_CPP_STATISTICS_H_ + +#include +#include +#include + +#include "ortools/math_opt/cpp/model.h" + +namespace operations_research::math_opt { + +// A range of values, first is the minimum, second is the maximum. +using Range = std::pair; + +// The ranges of the absolute values of the finite non-zero values in the model. +// +// Each range is optional since there may be no finite non-zero values +// (e.g. empty model, empty objective, all variables unbounded, ...). +struct ModelRanges { + // The linear and quadratic objective terms (not including the offset). + std::optional objective_terms; + + // The variables' lower and upper bounds. + std::optional variable_bounds; + + // The linear constraints' lower and upper bounds. + std::optional linear_constraint_bounds; + + // The coefficients of the variables in linear constraints. + std::optional linear_constraint_coefficients; +}; + +// Prints the ranges with std::setprecision(2) and std::scientific (restoring +// the stream flags after printing). +// +// It prints a multi-line table list of ranges. The last line does NOT end with +// a new line thus the caller should use std::endl if appropriate (see example). +// +// Example: +// +// std::cout << "Model xxx ranges:\n" << ComputeModelRanges(model) +// << std::endl; +// +// LOG(INFO) << "Model xxx ranges:\n" << ComputeModelRanges(model); +// +std::ostream& operator<<(std::ostream& out, const ModelRanges& ranges); + +// Returns the ranges of the finite non-zero values in the given model. +ModelRanges ComputeModelRanges(const Model& model); + +} // namespace operations_research::math_opt + +#endif // OR_TOOLS_MATH_OPT_CPP_STATISTICS_H_ diff --git a/ortools/math_opt/cpp/variable_and_expressions.h b/ortools/math_opt/cpp/variable_and_expressions.h index 8c3edee573..ace3dfb8a3 100644 --- a/ortools/math_opt/cpp/variable_and_expressions.h +++ b/ortools/math_opt/cpp/variable_and_expressions.h @@ -95,6 +95,7 @@ #include #include #include +#include #include #include "absl/container/flat_hash_map.h" @@ -231,8 +232,8 @@ class LinearExpression { // to this. // // Example: - // Variable a = ...; - // Variable b = ...; + // const Variable a = ...; + // const Variable b = ...; // const std::vector vars = {a, b}; // LinearExpression expr(8.0); // expr.AddSum(vars); @@ -251,6 +252,13 @@ class LinearExpression { template inline void AddSum(const Iterable& items); + // Creates a new LinearExpression object equal to the sum. The implementation + // is equivalent to: + // LinearExpression expr; + // expr.AddSum(items); + template + static inline LinearExpression Sum(const Iterable& items); + // Adds the inner product of left and right to this. // // Specifically, letting @@ -261,8 +269,8 @@ class LinearExpression { // to this. // // Example: - // Variable a = ...; - // Variable b = ...; + // const Variable a = ...; + // const Variable b = ...; // const std::vector left = {a, b}; // const std::vector right = {10.0, 2.0}; // LinearExpression expr(3.0); @@ -280,7 +288,7 @@ class LinearExpression { // Runtime requirements (or CHECK fails): // * left and right have an equal number of elements. // - // Note: The implementation is equivalent to: + // Note: The implementation is equivalent to the following pseudocode: // for(const auto& [l, r] : zip(left, right)) { // *this += l * r; // } @@ -291,6 +299,14 @@ class LinearExpression { inline void AddInnerProduct(const LeftIterable& left, const RightIterable& right); + // Creates a new LinearExpression object equal to the inner product. The + // implementation is equivalent to: + // LinearExpression expr; + // expr.AddInnerProduct(left, right); + template + static inline LinearExpression InnerProduct(const LeftIterable& left, + const RightIterable& right); + // Returns the terms in this expression. inline const VariableMap& terms() const; inline double offset() const; @@ -331,7 +347,7 @@ class LinearExpression { double offset_ = 0.0; }; -// Returns the sum of the elements of items. +// Returns the sum of the elements of items as a LinearExpression. // // Specifically, letting // (i_1, i_2, ..., i_n) = items @@ -339,8 +355,8 @@ class LinearExpression { // i_1 + i_2 + ... + i_n. // // Example: -// Variable a = ...; -// Variable b = ...; +// const Variable a = ...; +// const Variable b = ...; // const std::vector vars = {a, b, a}; // Sum(vars) // => 2.0 * a + b @@ -351,10 +367,13 @@ class LinearExpression { // expr.AddSum(items); // // See LinearExpression::AddSum() for a precise contract on the type Iterable. +// +// If the inner product cannot be represented as a LinearExpression, consider +// instead QuadraticExpression::Sum(). template inline LinearExpression Sum(const Iterable& items); -// Returns the inner product of left and right. +// Returns the inner product of left and right as a LinearExpression. // // Specifically, letting // (l_1, l_2 ..., l_n) = left, @@ -363,8 +382,8 @@ inline LinearExpression Sum(const Iterable& items); // l_1 * r_1 + l_2 * r_2 + ... + l_n * r_n. // // Example: -// Variable a = ...; -// Variable b = ...; +// const Variable a = ...; +// const Variable b = ...; // const std::vector left = {a, b}; // const std::vector right = {10.0, 2.0}; // InnerProduct(left, right); @@ -377,6 +396,9 @@ inline LinearExpression Sum(const Iterable& items); // // Requires that left and right have equal size, see // LinearExpression::AddInnerProduct for a precise contract on template types. +// +// If the inner product cannot be represented as a LinearExpression, consider +// instead QuadraticExpression::InnerProduct(). template inline LinearExpression InnerProduct(const LeftIterable& left, const RightIterable& right); @@ -628,6 +650,9 @@ class QuadraticTermKey { QuadraticProductId variable_ids_; }; +inline std::ostream& operator<<(std::ostream& ostr, + const QuadraticTermKey& key); + inline bool operator==(const QuadraticTermKey lhs, const QuadraticTermKey rhs); inline bool operator!=(const QuadraticTermKey lhs, const QuadraticTermKey rhs); @@ -738,6 +763,131 @@ class QuadraticExpression { inline QuadraticExpression& operator*=(double value); inline QuadraticExpression& operator/=(double value); + // Adds each element of items to this. + // + // Specifically, letting + // (i_1, i_2, ..., i_n) = items + // adds + // i_1 + i_2 + ... + i_n + // to this. + // + // Example: + // const Variable a = ...; + // const Variable b = ...; + // const std::vector vars = {a, b}; + // const std::vector terms = {2 * a * b}; + // QuadraticExpression expr = 8; + // expr.AddSum(vars); + // expr.AddSum(terms); + // Results in expr having the value 2 * a * b + a + b + 8.0. + // + // Compile time requirements: + // * Iterable is a sequence (an array or object with begin() and end()). + // * The type of an element of items is one of double, Variable, LinearTerm, + // LinearExpression, QuadraticTerm, or QuadraticExpression (or is + // implicitly convertible to one of these types, e.g. int). + // + // Note: The implementation is equivalent to: + // for(const auto item : items) { + // *this += item; + // } + template + inline void AddSum(const Iterable& items); + + // Returns the sum of the elements of items. + // + // Specifically, letting + // (i_1, i_2, ..., i_n) = items + // returns + // i_1 + i_2 + ... + i_n. + // + // Example: + // const Variable a = ...; + // const Variable b = ...; + // const std::vector terms = {a * a, 2 * a * b, 3 * b * a}; + // QuadraticExpression::Sum(vars) + // => a^2 + 5 a * b + // Note, instead of: + // QuadraticExpression expr(3.0); + // expr += QuadraticExpression::Sum(items); + // Prefer: + // expr.AddSum(items); + // + // See QuadraticExpression::AddSum() for a precise contract on the type + // Iterable. + template + static inline QuadraticExpression Sum(const Iterable& items); + + // Adds the inner product of left and right to this. + // + // Specifically, letting + // (l_1, l_2 ..., l_n) = left, + // (r_1, r_2, ..., r_n) = right, + // adds + // l_1 * r_1 + l_2 * r_2 + ... + l_n * r_n + // to this. + // + // Example: + // const Variable a = ...; + // const Variable b = ...; + // const std::vector vars = {a, b}; + // const std::vector coeffs = {10.0, 2.0}; + // QuadraticExpression expr = 3.0; + // expr.AddInnerProduct(coeffs, vars); + // expr.AddInnerProduct(vars, vars); + // Results in expr having the value a^2 + b^2 + 10.0 * a + 2.0 * b + 3.0. + // + // Compile time requirements: + // * LeftIterable and RightIterable are both sequences (arrays or objects + // with begin() and end()) + // * For both left and right, their elements are of type double, Variable, + // LinearTerm, LinearExpression, QuadraticTerm, or QuadraticExpression (or + // is implicitly convertible to one of these types, e.g. int). + // Runtime requirements (or CHECK fails): + // * The inner product value, and its constitutive intermediate terms, can be + // represented as a QuadraticExpression (potentially through an implicit + // conversion). + // * left and right have an equal number of elements. + // + // Note: The implementation is equivalent to the following pseudocode: + // for(const auto& [l, r] : zip(left, right)) { + // *this += l * r; + // } + // In particular, the multiplication will be performed on the types of the + // elements in left and right (take care with low precision types), but the + // addition will always use double precision. + template + inline void AddInnerProduct(const LeftIterable& left, + const RightIterable& right); + + // Returns the inner product of left and right. + // + // Specifically, letting + // (l_1, l_2 ..., l_n) = left, + // (r_1, r_2, ..., r_n) = right, + // returns + // l_1 * r_1 + l_2 * r_2 + ... + l_n * r_n. + // + // Example: + // const Variable a = ...; + // const Variable b = ...; + // const std::vector left = {a, a}; + // const std::vector left = {a, b}; + // QuadraticExpression::InnerProduct(left, right); + // -=> a^2 + a * b + // Note, instead of: + // QuadraticExpression expr(3.0); + // expr += QuadraticExpression::InnerProduct(left, right); + // Prefer: + // expr.AddInnerProduct(left, right); + // + // Requires that left and right have equal size, see + // QuadraticExpression::AddInnerProduct() for a precise contract on template + // types. + template + static inline QuadraticExpression InnerProduct(const LeftIterable& left, + const RightIterable& right); + // Compute the numeric value of this expression when variables are substituted // by their values in variable_values. // @@ -1217,23 +1367,30 @@ void LinearExpression::AddSum(const Iterable& items) { } template -LinearExpression Sum(const Iterable& items) { +LinearExpression LinearExpression::Sum(const Iterable& items) { LinearExpression result; result.AddSum(items); return result; } -template -void LinearExpression::AddInnerProduct(const LeftIterable& left, - const RightIterable& right) { +template +LinearExpression Sum(const Iterable& items) { + return LinearExpression::Sum(items); +} + +namespace internal { + +template +void AddInnerProduct(const LeftIterable& left, const RightIterable& right, + Expression& expr) { using std::begin; using std::end; auto l = begin(left); auto r = begin(right); - auto l_end = end(left); - auto r_end = end(right); + const auto l_end = end(left); + const auto r_end = end(right); for (; l != l_end && r != r_end; ++l, ++r) { - *this += (*l) * (*r); + expr += (*l) * (*r); } CHECK(l == l_end) << "left had more elements than right, sizes should be equal"; @@ -1241,14 +1398,28 @@ void LinearExpression::AddInnerProduct(const LeftIterable& left, << "right had more elements than left, sizes should be equal"; } +} // namespace internal + template -LinearExpression InnerProduct(const LeftIterable& left, - const RightIterable& right) { +void LinearExpression::AddInnerProduct(const LeftIterable& left, + const RightIterable& right) { + internal::AddInnerProduct(left, right, *this); +} + +template +LinearExpression LinearExpression::InnerProduct(const LeftIterable& left, + const RightIterable& right) { LinearExpression result; result.AddInnerProduct(left, right); return result; } +template +LinearExpression InnerProduct(const LeftIterable& left, + const RightIterable& right) { + return LinearExpression::InnerProduct(left, right); +} + const VariableMap& LinearExpression::terms() const { return terms_; } double LinearExpression::offset() const { return offset_; } @@ -1629,6 +1800,12 @@ H AbslHashValue(H h, const QuadraticTermKey& key) { key.typed_id().second.value(), key.storage()); } +std::ostream& operator<<(std::ostream& ostr, const QuadraticTermKey& key) { + ostr << "(" << Variable(key.storage(), key.typed_id().first) << ", " + << Variable(key.storage(), key.typed_id().second) << ")"; + return ostr; +} + bool operator==(const QuadraticTermKey lhs, const QuadraticTermKey rhs) { return lhs.storage() == rhs.storage() && lhs.typed_id() == rhs.typed_id(); } @@ -2229,6 +2406,34 @@ QuadraticExpression& QuadraticExpression::operator/=(const double value) { return *this; } +template +void QuadraticExpression::AddSum(const Iterable& items) { + for (const auto& item : items) { + *this += item; + } +} + +template +QuadraticExpression QuadraticExpression::Sum(const Iterable& items) { + QuadraticExpression result; + result.AddSum(items); + return result; +} + +template +void QuadraticExpression::AddInnerProduct(const LeftIterable& left, + const RightIterable& right) { + internal::AddInnerProduct(left, right, *this); +} + +template +QuadraticExpression QuadraticExpression::InnerProduct( + const LeftIterable& left, const RightIterable& right) { + QuadraticExpression result; + result.AddInnerProduct(left, right); + return result; +} + } // namespace math_opt } // namespace operations_research diff --git a/ortools/math_opt/io/BUILD.bazel b/ortools/math_opt/io/BUILD.bazel index 2d18725b74..76405497a3 100644 --- a/ortools/math_opt/io/BUILD.bazel +++ b/ortools/math_opt/io/BUILD.bazel @@ -38,3 +38,13 @@ cc_library( "@com_google_absl//absl/strings", ], ) + +cc_library( + name = "names_removal", + srcs = ["names_removal.cc"], + hdrs = ["names_removal.h"], + deps = [ + "//ortools/math_opt:model_cc_proto", + "//ortools/math_opt:model_update_cc_proto", + ], +) diff --git a/ortools/math_opt/io/names_removal.cc b/ortools/math_opt/io/names_removal.cc new file mode 100644 index 0000000000..b5c1bd020e --- /dev/null +++ b/ortools/math_opt/io/names_removal.cc @@ -0,0 +1,32 @@ +// Copyright 2010-2021 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "ortools/math_opt/io/names_removal.h" + +#include "ortools/math_opt/model.pb.h" +#include "ortools/math_opt/model_update.pb.h" + +namespace operations_research::math_opt { + +void RemoveNames(ModelProto& model) { + model.clear_name(); + model.mutable_variables()->clear_names(); + model.mutable_linear_constraints()->clear_names(); +} + +void RemoveNames(ModelUpdateProto& update) { + update.mutable_new_variables()->clear_names(); + update.mutable_new_linear_constraints()->clear_names(); +} + +} // namespace operations_research::math_opt diff --git a/ortools/math_opt/io/names_removal.h b/ortools/math_opt/io/names_removal.h new file mode 100644 index 0000000000..72fdbd134e --- /dev/null +++ b/ortools/math_opt/io/names_removal.h @@ -0,0 +1,36 @@ +// Copyright 2010-2021 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// Functions to remove names from models and updates. +// +// Input models can contain duplicated names which MathOpt solvers will +// refuse. The functions in this header can be use to mitigate that. +// +// These functions can also be used to anonymize models before saving them. +#ifndef OR_TOOLS_MATH_OPT_IO_NAMES_REMOVAL_H_ +#define OR_TOOLS_MATH_OPT_IO_NAMES_REMOVAL_H_ + +#include "ortools/math_opt/model.pb.h" +#include "ortools/math_opt/model_update.pb.h" + +namespace operations_research::math_opt { + +// Removes the model, variables and constraints names of the provided model. +void RemoveNames(ModelProto& model); + +// Removes the variables and constraints names of the provided update. +void RemoveNames(ModelUpdateProto& update); + +} // namespace operations_research::math_opt + +#endif // OR_TOOLS_MATH_OPT_IO_NAMES_REMOVAL_H_ diff --git a/ortools/math_opt/parameters.proto b/ortools/math_opt/parameters.proto index fa72e37fd8..ab8a4bfff5 100644 --- a/ortools/math_opt/parameters.proto +++ b/ortools/math_opt/parameters.proto @@ -115,6 +115,8 @@ enum EmphasisProto { } // Configures if potentially bad solver input is a warning or an error. +// +// TODO(b/196132970): implement this feature. message StrictnessProto { bool bad_parameter = 1; } diff --git a/ortools/math_opt/samples/basic_example.cc b/ortools/math_opt/samples/basic_example.cc index 34b0585873..e99a91615d 100644 --- a/ortools/math_opt/samples/basic_example.cc +++ b/ortools/math_opt/samples/basic_example.cc @@ -50,8 +50,8 @@ absl::Status Main() { switch (result.termination.reason) { case math_opt::TerminationReason::kOptimal: case math_opt::TerminationReason::kFeasible: - std::cout << "objective value: " << result.objective_value() << std::endl - << "value for variable x: " << result.variable_values().at(x) + std::cout << "Objective value: " << result.objective_value() << std::endl + << "Value for variable x: " << result.variable_values().at(x) << std::endl; return absl::OkStatus(); default: diff --git a/ortools/math_opt/samples/cocktail_hour.cc b/ortools/math_opt/samples/cocktail_hour.cc index ce524647e9..c7edabea58 100644 --- a/ortools/math_opt/samples/cocktail_hour.cc +++ b/ortools/math_opt/samples/cocktail_hour.cc @@ -28,7 +28,7 @@ // // In latex mode, the output can be piped directly to pdflatex, e.g. // blaze run -c opt \ -// ortools/math_opt/examples/cocktail_hour \ +// ortools/math_opt/examples/cpp/cocktail_hour \ // -- --num_ingredients 10 --mode latex | pdflatex -output-directory /tmp // will create a PDF in /tmp. #include diff --git a/ortools/math_opt/samples/facility_lp_benders.cc b/ortools/math_opt/samples/facility_lp_benders.cc index 49b63eb4ad..f5380376ce 100644 --- a/ortools/math_opt/samples/facility_lp_benders.cc +++ b/ortools/math_opt/samples/facility_lp_benders.cc @@ -658,12 +658,12 @@ absl::Status Main() { absl::Time start = absl::Now(); RETURN_IF_ERROR(FullProblem(instance, absl::GetFlag(FLAGS_solver_type))) << "full solve failed"; - std::cout << "Full solve time : " << absl::Now() - start << std::endl; + std::cout << "Full solve time: " << absl::Now() - start << std::endl; start = absl::Now(); RETURN_IF_ERROR(Benders(instance, absl::GetFlag(FLAGS_benders_precission), absl::GetFlag(FLAGS_solver_type))) << "Benders solve failed"; - std::cout << "Benders solve time : " << absl::Now() - start << std::endl; + std::cout << "Benders solve time: " << absl::Now() - start << std::endl; return absl::OkStatus(); } } // namespace diff --git a/ortools/math_opt/samples/integer_programming.cc b/ortools/math_opt/samples/integer_programming.cc index 2ced30b269..fccc22f452 100644 --- a/ortools/math_opt/samples/integer_programming.cc +++ b/ortools/math_opt/samples/integer_programming.cc @@ -58,7 +58,7 @@ absl::Status Main() { case math_opt::TerminationReason::kOptimal: case math_opt::TerminationReason::kFeasible: std::cout << "Problem solved in " << result.solve_time() << std::endl; - std::cout << "objective value: " << result.objective_value() << std::endl; + std::cout << "Objective value: " << result.objective_value() << std::endl; std::cout << "Variable values: [x=" << result.variable_values().at(x) << ", y=" << result.variable_values().at(y) << "]" << std::endl; return absl::OkStatus(); diff --git a/ortools/math_opt/samples/tsp.cc b/ortools/math_opt/samples/tsp.cc index b7757383f9..f7b6cd57c1 100644 --- a/ortools/math_opt/samples/tsp.cc +++ b/ortools/math_opt/samples/tsp.cc @@ -74,6 +74,8 @@ ABSL_FLAG(bool, test_instance, false, "Solve the test TSP instead of a random instance."); ABSL_FLAG(int, threads, 0, "How many threads to solve with, or solver default if <= 0."); +ABSL_FLAG(bool, solve_logs, false, + "Have the solver print logs to standard out."); namespace { @@ -260,6 +262,7 @@ absl::StatusOr SolveTsp( model.AddLinearConstraint(neighbors == 2, absl::StrCat("n_", i)); } math_opt::SolveArguments args; + args.parameters.enable_output = absl::GetFlag(FLAGS_solve_logs); const int threads = absl::GetFlag(FLAGS_threads); if (threads > 0) { args.parameters.threads = threads; diff --git a/ortools/math_opt/solvers/BUILD.bazel b/ortools/math_opt/solvers/BUILD.bazel index 05f9fe1b67..74b0d2716d 100644 --- a/ortools/math_opt/solvers/BUILD.bazel +++ b/ortools/math_opt/solvers/BUILD.bazel @@ -289,6 +289,7 @@ cc_library( "//ortools/math_opt/core:solver_interface", "//ortools/math_opt/core:sparse_submatrix", "//ortools/math_opt/core:sparse_vector_view", + "//ortools/math_opt/solvers/glpk:glpk_sparse_vector", "//ortools/math_opt/solvers/glpk:rays", "//ortools/math_opt/validators:callback_validator", "//ortools/port:proto_utils", diff --git a/ortools/math_opt/solvers/glpk/BUILD.bazel b/ortools/math_opt/solvers/glpk/BUILD.bazel index 0c5c0173e2..d10fd8d2bc 100644 --- a/ortools/math_opt/solvers/glpk/BUILD.bazel +++ b/ortools/math_opt/solvers/glpk/BUILD.bazel @@ -16,3 +16,10 @@ cc_library( "@glpk", ], ) + +cc_library( + name = "glpk_sparse_vector", + srcs = ["glpk_sparse_vector.cc"], + hdrs = ["glpk_sparse_vector.h"], + deps = ["//ortools/base"], +) diff --git a/ortools/math_opt/solvers/glpk/glpk_sparse_vector.cc b/ortools/math_opt/solvers/glpk/glpk_sparse_vector.cc new file mode 100644 index 0000000000..dafbd25350 --- /dev/null +++ b/ortools/math_opt/solvers/glpk/glpk_sparse_vector.cc @@ -0,0 +1,73 @@ +// Copyright 2010-2021 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "ortools/math_opt/solvers/glpk/glpk_sparse_vector.h" + +#include +#include + +#include "ortools/base/logging.h" + +namespace operations_research::math_opt { + +GlpkSparseVector::GlpkSparseVector(const int capacity) + : capacity_(capacity), + index_to_entry_(capacity + 1, kNotPresent), + indices_(capacity + 1, -1), + values_(capacity + 1, 0.0) { + CHECK_GE(capacity, 0); +} + +void GlpkSparseVector::Clear() { + // Resets the elements of the index_to_entry_ map we have modified. + for (int i = 1; i <= size_; ++i) { + index_to_entry_[indices_[i]] = kNotPresent; + } + + // Cleanup the used items to make sure we don't reuse those values by mistake + // later. + for (int i = 1; i <= size_; ++i) { + indices_[i] = -1; + values_[i] = 0.0; + } + + size_ = 0; +} + +void GlpkSparseVector::Load( + std::function getter) { + CHECK(getter != nullptr); + + Clear(); + + size_ = getter(indices_.data(), values_.data()); + + CHECK_GE(size_, 0); + CHECK_LE(size_, capacity_); + + // We don't know if the GLPK API has written to the first element but we reset + // those values anyway. + indices_[0] = -1; + values_[0] = 0.0; + + // Update index_to_entry_. + for (int i = 1; i <= size_; ++i) { + const int index = indices_[i]; + CHECK_GE(index, 1); + CHECK_LE(index, capacity_); + CHECK_EQ(index_to_entry_[index], kNotPresent) << "duplicated: " << index; + index_to_entry_[index] = i; + } +} + +} // namespace operations_research::math_opt diff --git a/ortools/math_opt/solvers/glpk/glpk_sparse_vector.h b/ortools/math_opt/solvers/glpk/glpk_sparse_vector.h new file mode 100644 index 0000000000..4c6161b39f --- /dev/null +++ b/ortools/math_opt/solvers/glpk/glpk_sparse_vector.h @@ -0,0 +1,201 @@ +// Copyright 2010-2021 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#ifndef OR_TOOLS_MATH_OPT_SOLVERS_GLPK_GLPK_SPARSE_VECTOR_H_ +#define OR_TOOLS_MATH_OPT_SOLVERS_GLPK_GLPK_SPARSE_VECTOR_H_ + +#include +#include +#include +#include + +#include "ortools/base/logging.h" + +namespace operations_research::math_opt { + +// Sparse vector in GLPK format. +// +// GLPK represents a sparse vector of size n with two arrays of size n+1, one +// for indices and one for values. The first element of each of these arrays is +// ignored (GLPK uses one-based indices). On top of that, the array of indices +// contains one-based indices (typically rows or columns indices). The entries +// are not necessarily sorted. +// +// For example to store a sparse vector where we have: +// +// idx | value +// ----+------ +// 1 | 2.5 +// 2 | +// 3 | -1.0 +// 4 | +// 5 | 0.5 +// +// GLPK would use two arrays: +// +// const int indices[] = { /*ignored*/-1, 3, 1, 5 }; +// const double values[] = { /*ignored*/NAN, -1.0, 2.5, 0.5 }; +// +// This class also keeps an additional vector which size is the capacity of the +// sparse vector (i.e. the corresponding size of a dense vector). It associates +// to each index an optional position of the corresponding entry in the indices +// and values arrays. This is used to make Set() and Get() O(1) and this makes +// Clear() O(size()) since indices associated to entries need to be cleared. +// +// This additional vector along with the ones used for indices and values are +// all pre-allocated to fit the capacity. Hence an instance of this class +// allocates: +// +// capacity * (2 * sizeof(int) + sizeof(double)) +// +// It is thus recommended to reuse the same instance multiple times instead of +// reallocating one for it to be efficient. +class GlpkSparseVector { + public: + // Builds a sparse vector with the provided capacity (i.e. the size of the + // vector it was dense). + // + // This operation has O(capacity) complexity (see the class documentation for + // allocated memory). + explicit GlpkSparseVector(const int capacity); + + // Returns the capacity (the size of the vector if it was dense). + int capacity() const { return capacity_; } + + // Returns the number of entries in the sparse vector. + int size() const { return size_; } + + // Returns the indices array of the GLPK sparse vector. + // + // Only values in [1, size()] are meaningful. + const int* indices() const { return indices_.data(); } + + // Returns the values array of the GLPK sparse vector. + // + // Only values in [1, size()] are meaningful. + const double* values() const { return values_.data(); } + + // Clears the sparse vector, removing all entries. + // + // This operation has O(size()) complexity. + void Clear(); + + // Returns the value at the given index if there is a corresponding entry or + // nullopt. + // + // It CHECKs that the index is in [1, capacity]. The operation has O(1) + // complexity. + inline std::optional Get(int index) const; + + // Changes the value of the given index, adding a new entry if necessary. + // + // Note that entries are only removed by Clear() or Load(). Setting a value to + // 0.0 does not remove the corresponding entry. + // + // It CHECKs that the index is in [1, capacity]. The operation has O(1) + // complexity. + inline void Set(int index, double value); + + // Replaces the content of the sparse vector by calling a GLPK API. + // + // Since GLPK functions have other parameters, here we expect the caller to + // provide a wrapping lambda expression that passes the indices and values + // buffers to the GLPK function and returns the number of written elements. + // + // It CHECKs that the returned number of elements is not greater than the + // capacity, that the indices are in [1, capacity] range and that there is no + // duplicated indices. + // + // Example: + // + // GlpkSparseVector row_values(num_cols); + // row_values.Set([&](int* const indices, double* const values) { + // return glp_get_mat_row(problem, row_index, indices, values); + // }); + // + void Load(std::function getter); + + private: + // Guard value used in index_to_entry_ to identify indices not in the sparse + // vector. + static constexpr int kNotPresent = std::numeric_limits::max(); + + // Capacity. + int capacity_; + + // Number of entries. + int size_ = 0; + + // For each dense index in [1, capacity], keeps the index of the corresponding + // entry in indices_ and values_. If the index i has a value in the sparse + // vector then indices_[index_to_entry_[i]] == i and + // values_[index_to_entry_[i]] is the corresponding value. If the index i does + // not have a value then index_to_entry_[i] == kNotPresent. + // + // Note that as for indices_ and values_, index_to_entry_[0] is unused. + std::vector index_to_entry_; + + // The GLPK one-based vector of entries' indices. Only values in [1, size_] + // are meaningful. + std::vector indices_; + + // The GLPK one-based vector of entries' values. Only values in [1, size_] are + // meaningful. + std::vector values_; +}; + +//////////////////////////////////////////////////////////////////////////////// +// Inline implementations +//////////////////////////////////////////////////////////////////////////////// + +std::optional GlpkSparseVector::Get(const int index) const { + CHECK_GE(index, 1); + CHECK_LE(index, capacity_); + + const int entry = index_to_entry_[index]; + if (entry == kNotPresent) { + return std::nullopt; + } + + DCHECK_GE(entry, 1); + DCHECK_LE(entry, capacity_); + DCHECK_EQ(indices_[entry], index); + + return values_[entry]; +} + +void GlpkSparseVector::Set(const int index, const double value) { + CHECK_GE(index, 1); + CHECK_LE(index, capacity_); + + const int entry = index_to_entry_[index]; + if (entry == kNotPresent) { + DCHECK_LT(size_, capacity_); + ++size_; + index_to_entry_[index] = size_; + indices_[size_] = index; + values_[size_] = value; + + return; + } + + DCHECK_GE(entry, 1); + DCHECK_LE(entry, capacity_); + DCHECK_EQ(indices_[entry], index); + + values_[entry] = value; +} + +} // namespace operations_research::math_opt + +#endif // OR_TOOLS_MATH_OPT_SOLVERS_GLPK_GLPK_SPARSE_VECTOR_H_ diff --git a/ortools/math_opt/solvers/glpk/rays.cc b/ortools/math_opt/solvers/glpk/rays.cc index d406c00938..b7a3645361 100644 --- a/ortools/math_opt/solvers/glpk/rays.cc +++ b/ortools/math_opt/solvers/glpk/rays.cc @@ -18,13 +18,13 @@ #include #include +#include "ortools/base/logging.h" #include "absl/status/status.h" #include "absl/status/statusor.h" #include "absl/strings/str_cat.h" -#include "ortools/base/logging.h" -#include "ortools/base/status_macros.h" #include "ortools/glpk/glpk_computational_form.h" #include "ortools/glpk/glpk_formatters.h" +#include "ortools/base/status_macros.h" namespace operations_research::math_opt { namespace { diff --git a/ortools/math_opt/solvers/glpk_solver.cc b/ortools/math_opt/solvers/glpk_solver.cc index 3cff3a26dc..c56ee7da04 100644 --- a/ortools/math_opt/solvers/glpk_solver.cc +++ b/ortools/math_opt/solvers/glpk_solver.cc @@ -57,6 +57,7 @@ #include "ortools/math_opt/parameters.pb.h" #include "ortools/math_opt/result.pb.h" #include "ortools/math_opt/solution.pb.h" +#include "ortools/math_opt/solvers/glpk/glpk_sparse_vector.h" #include "ortools/math_opt/solvers/glpk/rays.h" #include "ortools/math_opt/solvers/message_callback_data.h" #include "ortools/math_opt/sparse_containers.pb.h" @@ -833,18 +834,13 @@ double OffsetOnlyObjVal(glp_prob* const problem) { // glp_interior()). int OptStatus(glp_prob*) { return GLP_OPT; } -// Returns the error when a model or an update contains a quadratic objective. -absl::Status QuadraticObjectiveError() { - return absl::InvalidArgumentError( - "GLPK does not support quadratic objectives"); -} - } // namespace absl::StatusOr> GlpkSolver::New( const ModelProto& model, const InitArgs& init_args) { if (!model.objective().quadratic_coefficients().row_ids().empty()) { - return QuadraticObjectiveError(); + return absl::InvalidArgumentError( + "GLPK does not support quadratic objectives"); } return absl::WrapUnique(new GlpkSolver(model)); } @@ -1368,25 +1364,7 @@ void GlpkSolver::UpdateLinearConstraintMatrix( { // We reuse the same vectors for all calls to GLPK's API to limit // reallocations of these temporary buffers. - // - // We use constant size vectors to remove inefficiencies of - // std::vector::resize() that has linear cost in the size change (due to - // reset to 0 of values in the existing capacity on grow on resize). - // - // The value at index 0 is never used by GLPK's API (which are one-based). - std::vector data_indices(variables_.ids.size() + 1); - std::vector data_values(variables_.ids.size() + 1); - - // This shared vector (to prevent reallocation) will be used to store for - // each GLPK column the corresponding index in data_xxx vectors. The value - // kColNotPresent being used as a guard value to indicate that the column is - // not present. - // - // This uses the GLPK convention of ignoring the element 0 and using - // one-based indices. - constexpr int kColNotPresent = std::numeric_limits::max(); - std::vector col_to_data_element(variables_.ids.size() + 1, - kColNotPresent); + GlpkSparseVector data(static_cast(variables_.ids.size())); for (const auto& [row_id, row_coefficients] : SparseSubmatrixByRows(matrix_updates, /*start_row_id=*/0, @@ -1399,65 +1377,27 @@ void GlpkSolver::UpdateLinearConstraintMatrix( CHECK_EQ(linear_constraints_.ids[row_index - 1], row_id); // Read the current row coefficients. - const int initial_non_zeros = glp_get_mat_row( - problem_, row_index, data_indices.data(), data_values.data()); - CHECK_LE(initial_non_zeros + 1, data_indices.size()); - CHECK_LE(initial_non_zeros + 1, data_values.size()); - - // Update the col to data_xxx elements map. - for (int i = 1; i <= initial_non_zeros; ++i) { - col_to_data_element[data_indices[i]] = i; - } + data.Load([&](int* const indices, double* const values) { + return glp_get_mat_row(problem_, row_index, indices, values); + }); // Update the row data. - int non_zeros = initial_non_zeros; for (const auto [col_id, coefficient] : row_coefficients) { const int col_index = variables_.id_to_index.at(col_id); CHECK_EQ(variables_.ids[col_index - 1], col_id); - - // Here there are two cases: either a coefficient already exists for the - // given column, and we simply replace its value (potentially with a - // zero, which will remove the coefficient in the GLPK matrix), or we - // need to append it. - if (const int i = col_to_data_element[col_index]; i == kColNotPresent) { - ++non_zeros; - data_indices[non_zeros] = col_index; - data_values[non_zeros] = coefficient; - } else { - CHECK_EQ(data_indices[i], col_index); - data_values[i] = coefficient; - } + data.Set(col_index, coefficient); } - CHECK_LE(non_zeros + 1, data_indices.size()); - CHECK_LE(non_zeros + 1, data_values.size()); // Change the row values. - glp_set_mat_row(problem_, row_index, non_zeros, data_indices.data(), - data_values.data()); - - // Cleanup the used data_xxx items to make sure we don't reuse those - // values by mistake later. - for (int i = 1; i <= non_zeros; ++i) { - data_indices[i] = 0; - data_values[i] = 0; - } - - // Resets the elements of the map we modified. Here we ignore the new - // column indices that we have added in data_indices since those have - // not been put in the map. - for (int i = 1; i <= initial_non_zeros; ++i) { - col_to_data_element[data_indices[i]] = kColNotPresent; - } + glp_set_mat_row(problem_, row_index, data.size(), data.indices(), + data.values()); } } // Add new columns's coefficients of existing rows. The coefficients of new // columns in new rows will be added when adding new rows below. if (first_new_var_id.has_value()) { - // See the documentation for the existing rows above of the variables with - // the same names. - std::vector data_indices(linear_constraints_.ids.size() + 1); - std::vector data_values(linear_constraints_.ids.size() + 1); + GlpkSparseVector data(static_cast(linear_constraints_.ids.size())); for (const auto& [col_id, col_coefficients] : TransposeSparseSubmatrix( SparseSubmatrixByRows(matrix_updates, /*start_row_id=*/0, @@ -1471,37 +1411,22 @@ void GlpkSolver::UpdateLinearConstraintMatrix( // Prepare the column data replacing MathOpt ids by GLPK one-based row // indices. - int non_zeros = 0; + data.Clear(); for (const auto [row_id, coefficient] : MakeView(col_coefficients)) { const int row_index = linear_constraints_.id_to_index.at(row_id); CHECK_EQ(linear_constraints_.ids[row_index - 1], row_id); - - ++non_zeros; - data_indices[non_zeros] = row_index; - data_values[non_zeros] = coefficient; + data.Set(row_index, coefficient); } - CHECK_LE(non_zeros + 1, data_indices.size()); - CHECK_LE(non_zeros + 1, data_values.size()); // Change the column values. - glp_set_mat_col(problem_, col_index, non_zeros, data_indices.data(), - data_values.data()); - - // Cleanup the used data_xxx items to make sure we don't reuse those - // values by mistake later. - for (int i = 1; i <= non_zeros; ++i) { - data_indices[i] = 0; - data_values[i] = 0; - } + glp_set_mat_col(problem_, col_index, data.size(), data.indices(), + data.values()); } } // Add new rows, including the new columns' coefficients. if (first_new_cstr_id.has_value()) { - // See the documentation for the existing rows above of the variables with - // the same names. - std::vector data_indices(variables_.ids.size() + 1); - std::vector data_values(variables_.ids.size() + 1); + GlpkSparseVector data(static_cast(variables_.ids.size())); for (const auto& [row_id, row_coefficients] : SparseSubmatrixByRows(matrix_updates, /*start_row_id=*/*first_new_cstr_id, @@ -1515,28 +1440,16 @@ void GlpkSolver::UpdateLinearConstraintMatrix( // Prepare the row data replacing MathOpt ids by GLPK one-based column // indices. - int non_zeros = 0; + data.Clear(); for (const auto [col_id, coefficient] : row_coefficients) { const int col_index = variables_.id_to_index.at(col_id); CHECK_EQ(variables_.ids[col_index - 1], col_id); - - ++non_zeros; - data_indices[non_zeros] = col_index; - data_values[non_zeros] = coefficient; + data.Set(col_index, coefficient); } - CHECK_LE(non_zeros + 1, data_indices.size()); - CHECK_LE(non_zeros + 1, data_values.size()); // Change the row values. - glp_set_mat_row(problem_, row_index, non_zeros, data_indices.data(), - data_values.data()); - - // Cleanup the used data_xxx items to make sure we don't reuse those - // values by mistake later. - for (int i = 1; i <= non_zeros; ++i) { - data_indices[i] = 0; - data_values[i] = 0; - } + glp_set_mat_row(problem_, row_index, data.size(), data.indices(), + data.values()); } } } @@ -1649,13 +1562,6 @@ absl::Status GlpkSolver::AddPrimalOrDualRay( } absl::Status GlpkSolver::Update(const ModelUpdateProto& model_update) { - if (!model_update.objective_updates() - .quadratic_coefficients() - .row_ids() - .empty()) { - return QuadraticObjectiveError(); - } - // TODO(b/187027049): GLPK should not support modifying the model from another // thread (the allocation depends on the per-thread environment). We should // unit test that and see what is the actual behavior. If GLPK itself does not @@ -1735,11 +1641,10 @@ absl::Status GlpkSolver::Update(const ModelUpdateProto& model_update) { } bool GlpkSolver::CanUpdate(const ModelUpdateProto& model_update) { - // We return true even if we have a quadratic objective so that we don't force - // the caller to create a full model to get the error that quadratic - // objectives are not supported. The caller will get the correct error - // returned by GlpkSolver::Update(). - return true; + return model_update.objective_updates() + .quadratic_coefficients() + .row_ids() + .empty(); } MATH_OPT_REGISTER_SOLVER(SOLVER_TYPE_GLPK, GlpkSolver::New) diff --git a/ortools/math_opt/solvers/gurobi/g_gurobi.cc b/ortools/math_opt/solvers/gurobi/g_gurobi.cc index 77972159c2..a1070055dd 100644 --- a/ortools/math_opt/solvers/gurobi/g_gurobi.cc +++ b/ortools/math_opt/solvers/gurobi/g_gurobi.cc @@ -17,10 +17,10 @@ #include #include +#include "ortools/base/cleanup.h" #include "absl/status/status.h" #include "absl/status/statusor.h" #include "absl/strings/str_format.h" -#include "ortools/base/cleanup.h" #include "ortools/base/source_location.h" #include "ortools/base/status_builder.h" #include "ortools/base/status_macros.h" diff --git a/ortools/math_opt/solvers/gurobi/g_gurobi.h b/ortools/math_opt/solvers/gurobi/g_gurobi.h index 8007a76be3..7bc1941230 100644 --- a/ortools/math_opt/solvers/gurobi/g_gurobi.h +++ b/ortools/math_opt/solvers/gurobi/g_gurobi.h @@ -40,8 +40,9 @@ #include "absl/status/status.h" #include "absl/status/statusor.h" -#include "absl/types/span.h" #include "ortools/base/source_location.h" +#include "absl/types/span.h" + #include "ortools/gurobi/environment.h" namespace operations_research::math_opt { diff --git a/ortools/math_opt/solvers/gurobi_solver.cc b/ortools/math_opt/solvers/gurobi_solver.cc index bf47ab07b4..5a3b04d31a 100644 --- a/ortools/math_opt/solvers/gurobi_solver.cc +++ b/ortools/math_opt/solvers/gurobi_solver.cc @@ -751,12 +751,12 @@ absl::StatusOr GurobiSolver::ExtractSolveResultProto( const absl::Time start, const ModelSolveParametersProto& model_parameters) { SolveResultProto result; - // TODO(b/195295177): Add tests for rays in unbounded MIPs - RETURN_IF_ERROR(FillRays(model_parameters, result)); - ASSIGN_OR_RETURN((auto [solutions, solution_claims]), GetSolutions(model_parameters)); + // TODO(b/195295177): Add tests for rays in unbounded MIPs + RETURN_IF_ERROR(FillRays(model_parameters, solution_claims, result)); + for (auto& solution : solutions) { *result.add_solutions() = std::move(solution); } @@ -1147,10 +1147,14 @@ GurobiSolver::GetLpDualSolutionIfAvailable( absl::Status GurobiSolver::FillRays( const ModelSolveParametersProto& model_parameters, - SolveResultProto& result) { + const SolutionClaims solution_claims, SolveResultProto& result) { ASSIGN_OR_RETURN(const bool is_maximize, IsMaximize()); - if (gurobi_->IsAttrAvailable(GRB_DBL_ATTR_UNBDRAY) && - num_gurobi_variables_ > 0) { + // GRB_DBL_ATTR_UNBDRAY is sometimes incorrectly available for problems + // without variables. We also give priority to the conclusions obtained from + // dual solutions or bounds. + if (!solution_claims.dual_feasible_solution_exists && + num_gurobi_variables_ > 0 && + gurobi_->IsAttrAvailable(GRB_DBL_ATTR_UNBDRAY)) { ASSIGN_OR_RETURN(const std::vector grb_ray_var_values, gurobi_->GetDoubleAttrArray(GRB_DBL_ATTR_UNBDRAY, num_gurobi_variables_)); @@ -1159,8 +1163,12 @@ absl::Status GurobiSolver::FillRays( *primal_ray->mutable_variable_values(), model_parameters.variable_values_filter()); } - if (gurobi_->IsAttrAvailable(GRB_DBL_ATTR_FARKASDUAL) && - num_gurobi_constraints() + num_gurobi_variables_ > 0) { + // GRB_DBL_ATTR_FARKASDUAL is sometimes incorrectly available for problems + // without constraints. We also give priority to the conclusions obtained from + // primal solutions. + if (!solution_claims.primal_feasible_solution_exists && + num_gurobi_constraints() > 0 && + gurobi_->IsAttrAvailable(GRB_DBL_ATTR_FARKASDUAL)) { ASSIGN_OR_RETURN( DualRayProto dual_ray, GetGurobiDualRay(model_parameters.dual_values_filter(), @@ -1174,17 +1182,10 @@ absl::StatusOr GurobiSolver::GetQpSolution( const ModelSolveParametersProto& model_parameters) { ASSIGN_OR_RETURN((auto [primal_solution, found_primal_feasible_solution]), GetConvexPrimalSolutionIfAvailable(model_parameters)); - - // TODO(b/195295177): Update, seems GRB_DBL_ATTR_OBJBOUND is unavailable - // even for GRB_OPTIMAL, so the code below will only give a finite bound - // and a dual feasible status for GRB_OPTIMAL. - ASSIGN_OR_RETURN(const int grb_termination, - gurobi_->GetIntAttr(GRB_INT_ATTR_STATUS)); - bool dual_feasible_solution_exists = false; - ASSIGN_OR_RETURN(double best_dual_bound, GetBestDualBound()); - if (grb_termination == GRB_OPTIMAL || std::isfinite(best_dual_bound)) { - dual_feasible_solution_exists = true; - } + // TODO(b/225189115): Expand QpDualsTest to check maximization problems and + // other edge cases. + ASSIGN_OR_RETURN((auto [dual_solution, found_dual_feasible_solution]), + GetLpDualSolutionIfAvailable(model_parameters)); // Basis information is available when Gurobi uses QP simplex. As of v9.1 this // is not the default [1], so a user will need to explicitly set the Method // parameter in order for the following call to do anything interesting. @@ -1193,7 +1194,7 @@ absl::StatusOr GurobiSolver::GetQpSolution( const SolutionClaims solution_claims = { .primal_feasible_solution_exists = found_primal_feasible_solution, - .dual_feasible_solution_exists = dual_feasible_solution_exists}; + .dual_feasible_solution_exists = found_dual_feasible_solution}; if (!primal_solution.has_value() && !basis.has_value()) { return GurobiSolver::SolutionsAndClaims{.solution_claims = solution_claims}; @@ -1202,10 +1203,13 @@ absl::StatusOr GurobiSolver::GetQpSolution( SolutionProto& solution = solution_and_claims.solutions.emplace_back(SolutionProto()); if (primal_solution.has_value()) { - *solution.mutable_primal_solution() = std::move(*primal_solution); + *solution.mutable_primal_solution() = *std::move(primal_solution); + } + if (dual_solution.has_value()) { + *solution.mutable_dual_solution() = *std::move(dual_solution); } if (basis.has_value()) { - *solution.mutable_basis() = std::move(*basis); + *solution.mutable_basis() = *std::move(basis); } return solution_and_claims; } diff --git a/ortools/math_opt/solvers/gurobi_solver.h b/ortools/math_opt/solvers/gurobi_solver.h index 6e72773f6a..ae95d5d1f4 100644 --- a/ortools/math_opt/solvers/gurobi_solver.h +++ b/ortools/math_opt/solvers/gurobi_solver.h @@ -142,6 +142,7 @@ class GurobiSolver : public SolverInterface { absl::StatusOr ExtractSolveResultProto( absl::Time start, const ModelSolveParametersProto& model_parameters); absl::Status FillRays(const ModelSolveParametersProto& model_parameters, + const SolutionClaims solution_claims, SolveResultProto& result); absl::StatusOr GetSolutions( const ModelSolveParametersProto& model_parameters); diff --git a/ortools/math_opt/solvers/pdlp_solver.h b/ortools/math_opt/solvers/pdlp_solver.h index 8fb115e6f9..6b3ef4eb82 100644 --- a/ortools/math_opt/solvers/pdlp_solver.h +++ b/ortools/math_opt/solvers/pdlp_solver.h @@ -15,8 +15,6 @@ #define OR_TOOLS_MATH_OPT_SOLVERS_PDLP_SOLVER_H_ #include -#include -#include #include "absl/status/status.h" #include "absl/status/statusor.h" @@ -30,7 +28,6 @@ #include "ortools/math_opt/result.pb.h" #include "ortools/math_opt/solvers/pdlp_bridge.h" #include "ortools/pdlp/primal_dual_hybrid_gradient.h" -#include "ortools/pdlp/quadratic_program.h" #include "ortools/pdlp/solvers.pb.h" namespace operations_research { diff --git a/ortools/math_opt/tools/BUILD.bazel b/ortools/math_opt/tools/BUILD.bazel index 793973e2ba..2d2da997fb 100644 --- a/ortools/math_opt/tools/BUILD.bazel +++ b/ortools/math_opt/tools/BUILD.bazel @@ -9,7 +9,9 @@ cc_binary( "//ortools/math_opt:parameters_cc_proto", "//ortools/math_opt/core:solver_interface", "//ortools/math_opt/cpp:math_opt", + "//ortools/math_opt/cpp:statistics", "//ortools/math_opt/io:mps_converter", + "//ortools/math_opt/io:names_removal", "//ortools/math_opt/solvers:cp_sat_solver", "//ortools/math_opt/solvers:glop_solver", "//ortools/math_opt/solvers:glpk_solver", diff --git a/ortools/math_opt/tools/mathopt_solve_main.cc b/ortools/math_opt/tools/mathopt_solve_main.cc index fbfb422f8f..ee30e4f5a4 100644 --- a/ortools/math_opt/tools/mathopt_solve_main.cc +++ b/ortools/math_opt/tools/mathopt_solve_main.cc @@ -42,12 +42,18 @@ #include "ortools/base/status_macros.h" #include "ortools/math_opt/core/solver_interface.h" #include "ortools/math_opt/cpp/math_opt.h" +#include "ortools/math_opt/cpp/statistics.h" #include "ortools/math_opt/io/mps_converter.h" +#include "ortools/math_opt/io/names_removal.h" +#include "ortools/math_opt/io/proto_converter.h" #include "ortools/math_opt/parameters.pb.h" #include "ortools/util/status_macros.h" inline constexpr absl::string_view kMathOptBinaryFormat = "mathopt"; inline constexpr absl::string_view kMathOptTextFormat = "mathopt_txt"; +inline constexpr absl::string_view kLinearSolverBinaryFormat = "linear_solver"; +inline constexpr absl::string_view kLinearSolverTextFormat = + "linear_solver_txt"; inline constexpr absl::string_view kMPSFormat = "mps"; inline constexpr absl::string_view kAutoFormat = "auto"; @@ -73,17 +79,30 @@ struct SolverTypeProtoFormatter { ABSL_FLAG(std::string, input_file, "", "the file containing the model to solve; use --format to specify the " "file format"); -ABSL_FLAG(std::string, format, "auto", - absl::StrCat( - "the format of the --input_file; possible values:\n", "* ", - kMathOptBinaryFormat, ": for a MathOpt ModelProto in binary\n", - "* ", kMathOptTextFormat, ": when the proto is in text\n", "* ", - kMPSFormat, ": for MPS file (which can be GZiped)\n", "* ", - kAutoFormat, ": to guess the format from the file extension:\n", - " - '", kPbExt, "', '", kProtoExt, "': ", kMathOptBinaryFormat, - "\n", " - '", kPbTxtExt, "', '", kTextProtoExt, - "': ", kMathOptTextFormat, "\n", " - '", kMPSExt, "', '", - kMPSGzipExt, "': ", kMPSFormat)); +ABSL_FLAG( + std::string, format, "auto", + absl::StrCat( + "the format of the --input_file; possible values:\n", + // + "* ", kMathOptBinaryFormat, ": for a MathOpt ModelProto in binary\n", + // + "* ", kMathOptTextFormat, ": when the proto is in text\n", + // + "* ", kLinearSolverBinaryFormat, + ": for a LinearSolver MPModelProto in binary\n", + // + "* ", kLinearSolverTextFormat, ": when the proto is in text\n", + // + "* ", kMPSFormat, ": for MPS file (which can be GZiped)\n", + // + "* ", kAutoFormat, ": to guess the format from the file extension:\n", + // + " - '", kPbExt, "', '", kProtoExt, "': ", kMathOptBinaryFormat, "\n", + // + " - '", kPbTxtExt, "', '", kTextProtoExt, "': ", kMathOptTextFormat, + "\n", + // + " - '", kMPSExt, "', '", kMPSGzipExt, "': ", kMPSFormat)); ABSL_FLAG( std::vector, update_files, {}, absl::StrCat( @@ -102,6 +121,12 @@ ABSL_FLAG(bool, solver_logs, false, "use a message callback to print the solver convergence logs"); ABSL_FLAG(absl::Duration, time_limit, absl::InfiniteDuration(), "the time limit to use for the solve"); +ABSL_FLAG(bool, names, true, + "use the names in the input models; ignoring names is useful when " + "the input contains duplicates"); +ABSL_FLAG(bool, ranges, false, + "prints statistics about the ranges of the model values"); +ABSL_FLAG(bool, print_model, false, "prints the model to stdout"); namespace operations_research { namespace math_opt { @@ -138,6 +163,15 @@ absl::StatusOr ReadModel(const absl::string_view file_path, if (format == kMathOptTextFormat) { return file::GetTextProto(file_path, file::Defaults()); } + if (format == kLinearSolverBinaryFormat || + format == kLinearSolverTextFormat) { + ASSIGN_OR_RETURN( + MPModelProto linear_solver_model, + format == kLinearSolverBinaryFormat + ? file::GetBinaryProto(file_path, file::Defaults()) + : file::GetTextProto(file_path, file::Defaults())); + return MPModelProtoToMathOptModel(linear_solver_model); + } if (format == kMPSFormat) { return ReadMpsFile(file_path); } @@ -159,6 +193,28 @@ absl::StatusOr ReadModelUpdate( absl::StrCat("invalid format in ReadModelUpdate(): ", format)); } +// Prints the model. +void PrintModel(Model& model) { + if (model.is_maximize()) { + std::cout << "max"; + } else { + std::cout << "min"; + } + std::cout << ' ' << model.ObjectiveAsQuadraticExpression() << std::endl; + + std::cout << "variables:" << std::endl; + for (const Variable v : model.SortedVariables()) { + std::cout << ' ' << v.lower_bound() << " ≤ " << v << " ≤ " + << v.upper_bound() << std::endl; + } + + std::cout << "constraints:" << std::endl; + for (const LinearConstraint c : model.SortedLinearConstraints()) { + std::cout << ' ' << c << ": " << model.AsBoundedLinearExpression(c) + << std::endl; + } +} + // Prints the summary of the solve result. absl::Status PrintSummary(const SolveResult& result) { std::cout << "Solve finished:\n" @@ -212,7 +268,7 @@ absl::Status RunSolver() { << "."; } - OR_ASSIGN_OR_RETURN3(const ModelProto model_proto, + OR_ASSIGN_OR_RETURN3(ModelProto model_proto, ReadModel(input_file_path, format), _ << "failed to read " << input_file_path); @@ -223,7 +279,14 @@ absl::Status RunSolver() { model_updates.emplace_back(std::move(update)); } - // Solve the problem. + if (!absl::GetFlag(FLAGS_names)) { + RemoveNames(model_proto); + for (ModelUpdateProto& update : model_updates) { + RemoveNames(update); + } + } + + // Parse the problem and the updates. ASSIGN_OR_RETURN(const std::unique_ptr model, Model::FromModelProto(model_proto)); for (int u = 0; u < model_updates.size(); ++u) { @@ -232,6 +295,17 @@ absl::Status RunSolver() { << "failed to apply the update file: " << update_file_paths[u]; } + if (absl::GetFlag(FLAGS_ranges)) { + std::cout << "Ranges of finite non-zero values in the model:\n" + << ComputeModelRanges(*model) << std::endl; + } + + // Optionally prints the problem. + if (absl::GetFlag(FLAGS_print_model)) { + PrintModel(*model); + } + + // Solve the problem. SolveArguments solve_args = { .parameters = {.time_limit = absl::GetFlag(FLAGS_time_limit)}, }; diff --git a/ortools/math_opt/validators/result_validator.cc b/ortools/math_opt/validators/result_validator.cc index b88bf463bc..c404c5ed4c 100644 --- a/ortools/math_opt/validators/result_validator.cc +++ b/ortools/math_opt/validators/result_validator.cc @@ -181,7 +181,6 @@ absl::Status CheckDualSolutionAndStatusConsistency( // Assumes ValidateTermination has been called and ValidateProblemStatusProto // has been called on result.solve_stats.problem_status. -// TODO(b/212685946): add ray checks absl::Status ValidateTerminationConsistency(const SolveResultProto& result) { const ProblemStatusProto status = result.solve_stats().problem_status(); switch (result.termination().reason()) { @@ -286,12 +285,26 @@ absl::Status ValidateResult(const SolveResultProto& result, RETURN_IF_ERROR( ValidateSolutions(result.solutions(), parameters, model_summary)); + if (result.primal_rays_size() > 0 && + result.solve_stats().problem_status().dual_status() == + FEASIBILITY_STATUS_FEASIBLE) { + return absl::InvalidArgumentError( + "solve_stats.problem_status.dual_status = FEASIBILITY_STATUS_FEASIBLE, " + "but a primal ray is returned"); + } for (int i = 0; i < result.primal_rays_size(); ++i) { RETURN_IF_ERROR(ValidatePrimalRay(result.primal_rays(i), parameters.variable_values_filter(), model_summary)) << "Invalid primal_rays[" << i << "]"; } + if (result.dual_rays_size() > 0 && + result.solve_stats().problem_status().primal_status() == + FEASIBILITY_STATUS_FEASIBLE) { + return absl::InvalidArgumentError( + "solve_stats.problem_status.primal_status = " + "FEASIBILITY_STATUS_FEASIBLE, but a dual ray is returned"); + } for (int i = 0; i < result.dual_rays_size(); ++i) { RETURN_IF_ERROR( ValidateDualRay(result.dual_rays(i), parameters, model_summary))