From 2f63f6540f903b56d55d10f4dcfc675da50fe79d Mon Sep 17 00:00:00 2001 From: Ulfw Date: Mon, 7 Oct 2024 21:44:48 +0200 Subject: [PATCH 01/10] Add MOSEK support for linear interface and math_opt --- CMakeLists.txt | 5 + cmake/FindMOSEK.cmake | 171 ++ cmake/check_deps.cmake | 5 + cmake/cpp.cmake | 4 + cmake/system_deps.cmake | 4 + examples/cpp/integer_programming.cc | 1 + examples/cpp/linear_programming.cc | 1 + ...trawberry_fields_with_column_generation.cc | 6 + .../cpp/uncapacitated_facility_location.cc | 5 + ortools/linear_solver/BUILD.bazel | 19 + ortools/linear_solver/CMakeLists.txt | 4 +- ortools/linear_solver/linear_solver.cc | 19 + ortools/linear_solver/linear_solver.h | 7 + ortools/linear_solver/linear_solver.proto | 3 + ortools/linear_solver/mosek_interface.cc | 1402 +++++++++++++++++ ortools/math_opt/cpp/parameters.h | 2 + ortools/math_opt/parameters.proto | 6 + ortools/math_opt/python/parameters.py | 3 + ortools/math_opt/solvers/BUILD.bazel | 53 + ortools/math_opt/solvers/CMakeLists.txt | 6 + ortools/math_opt/solvers/mosek.proto | 28 + ortools/math_opt/solvers/mosek/BUILD.bazel | 39 + ortools/math_opt/solvers/mosek/mosekwrp.cc | 546 +++++++ ortools/math_opt/solvers/mosek/mosekwrp.h | 180 +++ ortools/math_opt/solvers/mosek_solver.cc | 1181 ++++++++++++++ ortools/math_opt/solvers/mosek_solver.h | 188 +++ ortools/service/v1/mathopt/parameters.proto | 5 + 27 files changed, 3892 insertions(+), 1 deletion(-) create mode 100644 cmake/FindMOSEK.cmake create mode 100644 ortools/linear_solver/mosek_interface.cc create mode 100644 ortools/math_opt/solvers/mosek.proto create mode 100644 ortools/math_opt/solvers/mosek/BUILD.bazel create mode 100644 ortools/math_opt/solvers/mosek/mosekwrp.cc create mode 100644 ortools/math_opt/solvers/mosek/mosekwrp.h create mode 100644 ortools/math_opt/solvers/mosek_solver.cc create mode 100644 ortools/math_opt/solvers/mosek_solver.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 0fd2682a29b..e59aae723fc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -268,6 +268,11 @@ endif() CMAKE_DEPENDENT_OPTION(USE_GUROBI "Use the Gurobi solver" ON "BUILD_CXX" OFF) message(STATUS "Gurobi support: ${USE_GUROBI}") +## MOSEK +# see: https://mosek.com +CMAKE_DEPENDENT_OPTION(USE_HIGHS "Use the MOSEK solver" ON "BUILD_CXX" OFF) +message(STATUS "MOSEK support: ${USE_MOSEK}") + ## HiGHS # see: https://github.com/ERGO-Code/HiGHS CMAKE_DEPENDENT_OPTION(USE_HIGHS "Use the HiGHS solver" ON "BUILD_CXX" OFF) diff --git a/cmake/FindMOSEK.cmake b/cmake/FindMOSEK.cmake new file mode 100644 index 00000000000..515c292266a --- /dev/null +++ b/cmake/FindMOSEK.cmake @@ -0,0 +1,171 @@ +# Copyright 2010-2024 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. + +#[=======================================================================[.rst: +FindMosek +-------- + +This module determines the Mosek library of the system. + +IMPORTED Targets +^^^^^^^^^^^^^^^^ + +This module defines :prop_tgt:`IMPORTED` target ``mosek::mosek``, if +Mosek has been found. + +Result Variables +^^^^^^^^^^^^^^^^ + +This module defines the following variables: + +:: + +MOSEK_FOUND - True if Mosek found. + +Hints +^^^^^ + +A user may set ``MOSEK_PLATFORM_DIR`` to a Mosek installation platform +directoru to tell this module where to look, or ``MOSEK_BASE`` to point +to path where the ``mosek/`` directory is located. +#]=======================================================================] +set(MOSEK_FOUND FALSE) +message(STATUS "Locating MOSEK") + +if(CMAKE_C_COMPILER_LOADED) + include (CheckIncludeFile) + include (CheckCSourceCompiles) +elseif(CMAKE_CXX_COMPILER_LOADED) + include (CheckIncludeFileCXX) + include (CheckCXXSourceCompiles) +else() + message(FATAL_ERROR "FindMosek only works if either C or CXX language is enabled") +endif() + +if(APPLE) + if(CMAKE_SYSTEM_PROCESSOR MATCHES "^(aarch64|arm64)") + SET(MOSEK_PLATFORM_NAME "osxaarch64") + elseif() + message(FATAL_ERROR "Mosek not supported for ${CMAKE_SYSTEM} / ${CMAKE_SYSTEM_PROCESSOR}") + endif() +elseif(UNIX) + if(CMAKE_SYSTEM_PROCESSOR MATCHES "^(aarch64|arm64)") + SET(MOSEK_PLATFORM_NAME "linuxaarch64") + else() + SET(MOSEK_PLATFORM_NAME "linux64x86") + endif() +elseif(MSVC) + SET(MOSEK_PLATFORM_NAME "win64x86") +else() + message(FATAL_ERROR "Mosek not supported for ${CMAKE_SYSTEM}") +endif() + +function(FindMosekPlatformInPath RESULT PATH) + SET(${RESULT} "") + if(EXISTS "${PATH}/mosek") + SET(dirlist "") + FILE(GLOB entries LIST_DIRECTORIES true "${PATH}/mosek/*") + FOREACH(f ${entries}) + if(IS_DIRECTORY "${f}") + get_filename_component(bn "${f}" NAME) + if("${bn}" MATCHES "^[0-9]+[.][0-9]+$") + if (${bn} GREATER_EQUAL "10.0") + LIST(APPEND dirlist "${bn}") + endif() + endif() + endif() + ENDFOREACH() + LIST(SORT dirlist COMPARE NATURAL ORDER DESCENDING) + + if(MOSEK_PLATFORM_NAME) + foreach(MOSEK_VERSION ${dirlist}) + SET(MSKPFDIR "${PATH}/mosek/${MOSEK_VERSION}/tools/platform/${MOSEK_PLATFORM_NAME}") + if(EXISTS "${MSKPFDIR}") + SET(${RESULT} ${MSKPFDIR} PARENT_SCOPE) + return() + endif() + endforeach() + endif() + endif() +endfunction() + + +# Where to look for MOSEK: +# Standard procedure in Linux/OSX is to install MOSEK in the home directory, i.e. +# $HOME/mosek/X.Y/... +# Option 1. The user can specify when running CMake where the MOSEK platform directory is located, e.g. +# -DMOSEK_PLATFORM_DIR=$HOME/mosek/10.2/tools/platform/linux64x86/ +# in which case no search is performed. +# Option 2. The user can specify MOSEK_ROOT when running cmake. MOSEK_ROOT is +# the directory where the root mosek/ tree is located. +# Option 3. Automatic search. We will then attempt to search in the default +# locations, and if that fails, assume it is installed in a system location. +# For option 2 and 3, the newest MOSEK version will be chosen if more are available. + +if(MOSEK_PLATFORM_DIR) + # User defined platform dir directly +elseif(MOSEK_BASE) + # Look under for MOSEK_ROOT/X.Y + FindMosekPlatformInPath(MOSEK_PLATFORM_DIR "${MOSEK_BASE}") + if(NOT MOSEK_PLATFORM_DIR) + message(FATAL_ERROR " Could not locate MOSEK platform directory under ${MOSEK_BASE}") + endif() +endif() +if(NOT MOSEK_PLATFORM_DIR) + # Look in users home dir + if(EXISTS $ENV{HOME}) + FindMosekPlatformInPath(MOSEK_PLATFORM_DIR "$ENV{HOME}") + endif() +endif() +if(NOT MOSEK_PLATFORM_DIR) + # Look in users home dir + if(ENV{HOMEDRIVE} AND ENV{HOMEPATH}) + FindMosekPlatformInPath(MOSEK_PLATFORM_DIR"$ENV{HOMEDRIVE}$ENV{HOMEPATH}") + endif() +endif() + +if(NOT MOSEK_PLATFORM_DIR) + message(FATAL_ERROR " MOSEK_PLATFORM_DIR could not be detected") +else() + message(STATUS " MOSEK_PLATFORM_DIR detected: ${MOSEK_PLATFORM_DIR}") + set(MOSEK_PFDIR_FOUND TRUE) +endif() + + +if(MOSEK_PFDIR_FOUND AND NOT TARGET Mosek::Mosek) + add_library(mosek::mosek UNKNOWN IMPORTED) + + find_path(MOSEKINC mosek.h HINTS ${MOSEK_PLATFORM_DIR}/h) + find_library(LIBMOSEK mosek64 HINTS ${MOSEK_PLATFORM_DIR}/bin) + + if(LIBMOSEK) + set_target_properties(mosek::mosek PROPERTIES IMPORTED_LOCATION ${LIBMOSEK}) + endif() + + if (MOSEKINC) + target_include_directories(mosek::mosek INTERFACE "${MOSEKINC}") + endif() +elseif(NOT TARGET Mosek::Mosek) + add_library(mosek::mosek UNKNOWN IMPORTED) + + find_path(MOSEKINC mosek.h) + find_library(LIBMOSEK mosek64) + + if(LIBMOSEK) + set_target_properties(mosek::mosek PROPERTIES IMPORTED_LOCATION ${LIBMOSEK}) + endif() + + if (MOSEKINC) + target_include_directories(mosek::mosek INTERFACE "${MOSEKINC}") + endif() +endif() diff --git a/cmake/check_deps.cmake b/cmake/check_deps.cmake index 88bc8c4db24..a0b4ba90b86 100644 --- a/cmake/check_deps.cmake +++ b/cmake/check_deps.cmake @@ -116,6 +116,11 @@ if(USE_CPLEX AND NOT TARGET CPLEX::CPLEX) message(FATAL_ERROR "Target CPLEX::CPLEX not available.") endif() +# Check optional Dependencies +if(USE_MOSEK AND NOT TARGET mosek::mosek) + message(FATAL_ERROR "Target mosek::mosek not available.") +endif() + # CXX Test if(BUILD_TESTING) if(NOT TARGET GTest::gtest_main) diff --git a/cmake/cpp.cmake b/cmake/cpp.cmake index 6f5788e8510..692de2bfb80 100644 --- a/cmake/cpp.cmake +++ b/cmake/cpp.cmake @@ -90,6 +90,9 @@ endif() if(USE_CPLEX) list(APPEND OR_TOOLS_COMPILE_DEFINITIONS "USE_CPLEX") endif() +if(USE_MOSEK) + list(APPEND OR_TOOLS_COMPILE_DEFINITIONS "USE_MOSEK") +endif() if(WIN32) list(APPEND OR_TOOLS_COMPILE_DEFINITIONS "__WIN32__") @@ -565,6 +568,7 @@ target_link_libraries(${PROJECT_NAME} PUBLIC ${CPLEX_DEPS} ${GLPK_DEPS} ${HIGHS_DEPS} + ${MOSEK_DEPS} ${PDLP_DEPS} ${SCIP_DEPS} Threads::Threads) diff --git a/cmake/system_deps.cmake b/cmake/system_deps.cmake index 4dd0402b6fa..b27c5f2023a 100644 --- a/cmake/system_deps.cmake +++ b/cmake/system_deps.cmake @@ -87,6 +87,10 @@ if(USE_CPLEX) find_package(CPLEX REQUIRED) endif() +if(USE_MOSEK) + find_package(MOSEK REQUIRED) +endif() + # CXX Test if(BUILD_TESTING AND NOT BUILD_googletest) find_package(GTest REQUIRED) diff --git a/examples/cpp/integer_programming.cc b/examples/cpp/integer_programming.cc index 33e98eae682..c13fcb13291 100644 --- a/examples/cpp/integer_programming.cc +++ b/examples/cpp/integer_programming.cc @@ -88,6 +88,7 @@ void RunAllExamples() { RunIntegerProgrammingExample("GLPK"); RunIntegerProgrammingExample("CPLEX"); RunIntegerProgrammingExample("XPRESS"); + RunIntegerProgrammingExample("MOSEK"); } } // namespace operations_research diff --git a/examples/cpp/linear_programming.cc b/examples/cpp/linear_programming.cc index 94639e3483e..05c2dcee811 100644 --- a/examples/cpp/linear_programming.cc +++ b/examples/cpp/linear_programming.cc @@ -118,6 +118,7 @@ void RunAllExamples() { RunLinearProgrammingExample("XPRESS_LP"); RunLinearProgrammingExample("PDLP"); RunLinearProgrammingExample("HIGHS"); + RunLinearProgrammingExample("MOSEK_LP"); } } // namespace operations_research diff --git a/examples/cpp/strawberry_fields_with_column_generation.cc b/examples/cpp/strawberry_fields_with_column_generation.cc index c9f376b1958..eae4c715b30 100644 --- a/examples/cpp/strawberry_fields_with_column_generation.cc +++ b/examples/cpp/strawberry_fields_with_column_generation.cc @@ -633,6 +633,12 @@ int main(int argc, char** argv) { solver_type = operations_research::MPSolver::CPLEX_LINEAR_PROGRAMMING; found = true; } +#endif +#if defined(USE_MOSEK) + if (absl::GetFlag(FLAGS_colgen_solver) == "mosek") { + solver_type = operations_research::MPSolver::MOSEK_LINEAR_PROGRAMMING; + found = true; + } #endif if (!found) { LOG(ERROR) << "Unknown solver " << absl::GetFlag(FLAGS_colgen_solver); diff --git a/examples/cpp/uncapacitated_facility_location.cc b/examples/cpp/uncapacitated_facility_location.cc index e91707e0c45..3514928d624 100644 --- a/examples/cpp/uncapacitated_facility_location.cc +++ b/examples/cpp/uncapacitated_facility_location.cc @@ -219,6 +219,11 @@ void RunAllExamples(int32_t facilities, int32_t clients, double fix_cost) { LOG(INFO) << "---- Integer programming example with CPLEX ----"; UncapacitatedFacilityLocation(facilities, clients, fix_cost, MPSolver::CPLEX_MIXED_INTEGER_PROGRAMMING); +#endif // USE_CPLEX +#if defined(USE_MOSEK) + LOG(INFO) << "---- Integer programming example with MOSEK ----"; + UncapacitatedFacilityLocation(facilities, clients, fix_cost, + MPSolver::MOSEK_MIXED_INTEGER_PROGRAMMING); #endif // USE_CPLEX LOG(INFO) << "---- Integer programming example with CP-SAT ----"; UncapacitatedFacilityLocation(facilities, clients, fix_cost, diff --git a/ortools/linear_solver/BUILD.bazel b/ortools/linear_solver/BUILD.bazel index 618e1921a6c..5276e51c442 100644 --- a/ortools/linear_solver/BUILD.bazel +++ b/ortools/linear_solver/BUILD.bazel @@ -143,6 +143,19 @@ config_setting( }, ) +bool_flag( + name = "with_mosek", + build_setting_default = False, +) + +config_setting( + name = "use_mosek", + flag_values = { + ":with_mosek": "true", + }, +) + + proto_library( name = "linear_solver_proto", srcs = ["linear_solver.proto"], @@ -209,6 +222,9 @@ cc_library( }) + select({ ":use_cplex": ["cplex_interface.cc"], "//conditions:default": [], + }) + select({ + ":use_mosek": ["mosek_interface.cc"], + "//conditions:default": [], }), hdrs = [ "linear_expr.h", @@ -251,6 +267,9 @@ cc_library( }) + select({ ":use_cplex": ["-DUSE_CPLEX"], "//conditions:default": [], + }) + select({ + ":use_mosek": ["-DUSE_MOSEK"], + "//conditions:default": [], }), deps = [ ":linear_solver_cc_proto", diff --git a/ortools/linear_solver/CMakeLists.txt b/ortools/linear_solver/CMakeLists.txt index 74d20df48bc..cfbe3e1469d 100644 --- a/ortools/linear_solver/CMakeLists.txt +++ b/ortools/linear_solver/CMakeLists.txt @@ -31,7 +31,8 @@ set_target_properties(${NAME} PROPERTIES ) target_include_directories(${NAME} PUBLIC $ - $) + $ + $<$:${MOSEK_PLATFORM_DIR}/h>) target_link_libraries(${NAME} PRIVATE absl::memory absl::strings @@ -41,6 +42,7 @@ target_link_libraries(${NAME} PRIVATE $<$:Coin::Cbc> $<$:Coin::Clp> $<$:CPLEX::CPLEX> + $<$:mosek::mosek> $<$:GLPK::GLPK> $<$:highs::highs> $<$:Eigen3::Eigen> diff --git a/ortools/linear_solver/linear_solver.cc b/ortools/linear_solver/linear_solver.cc index 77de096ce99..9a93fcf3547 100644 --- a/ortools/linear_solver/linear_solver.cc +++ b/ortools/linear_solver/linear_solver.cc @@ -92,6 +92,7 @@ bool SolverTypeIsMip(MPModelRequest::SolverType solver_type) { case MPModelRequest::HIGHS_LINEAR_PROGRAMMING: case MPModelRequest::XPRESS_LINEAR_PROGRAMMING: case MPModelRequest::CPLEX_LINEAR_PROGRAMMING: + case MPModelRequest::MOSEK_LINEAR_PROGRAMMING: return false; case MPModelRequest::SCIP_MIXED_INTEGER_PROGRAMMING: @@ -104,6 +105,7 @@ bool SolverTypeIsMip(MPModelRequest::SolverType solver_type) { case MPModelRequest::HIGHS_MIXED_INTEGER_PROGRAMMING: case MPModelRequest::XPRESS_MIXED_INTEGER_PROGRAMMING: case MPModelRequest::CPLEX_MIXED_INTEGER_PROGRAMMING: + case MPModelRequest::MOSEK_MIXED_INTEGER_PROGRAMMING: return true; } LOG(DFATAL) << "Invalid SolverType: " << solver_type; @@ -404,6 +406,9 @@ extern MPSolverInterface* BuildGurobiInterface(bool mip, #if defined(USE_CPLEX) extern MPSolverInterface* BuildCplexInterface(bool mip, MPSolver* const solver); #endif +#if defined(USE_MOSEK) +extern MPSolverInterface* BuildMosekInterface(bool mip, MPSolver* const solver); +#endif extern MPSolverInterface* BuildXpressInterface(bool mip, MPSolver* const solver); @@ -458,6 +463,12 @@ MPSolverInterface* BuildSolverInterface(MPSolver* const solver) { return BuildCplexInterface(false, solver); case MPSolver::CPLEX_MIXED_INTEGER_PROGRAMMING: return BuildCplexInterface(true, solver); +#endif +#if defined(USE_MOSEK) + case MPSolver::MOSEK_LINEAR_PROGRAMMING: + return BuildMosekInterface(false, solver); + case MPSolver::MOSEK_MIXED_INTEGER_PROGRAMMING: + return BuildMosekInterface(true, solver); #endif case MPSolver::XPRESS_MIXED_INTEGER_PROGRAMMING: return BuildXpressInterface(true, solver); @@ -542,6 +553,12 @@ bool MPSolver::SupportsProblemType(OptimizationProblemType problem_type) { problem_type == CPLEX_MIXED_INTEGER_PROGRAMMING) { return true; } +#endif +#ifdef USE_MOSEK + if (problem_type == MOSEK_LINEAR_PROGRAMMING || + problem_type == MOSEK_MIXED_INTEGER_PROGRAMMING) { + return true; + } #endif if (problem_type == XPRESS_MIXED_INTEGER_PROGRAMMING || problem_type == XPRESS_LINEAR_PROGRAMMING) { @@ -573,6 +590,7 @@ constexpr {MPSolver::GLPK_LINEAR_PROGRAMMING, "glpk_lp"}, {MPSolver::HIGHS_LINEAR_PROGRAMMING, "highs_lp"}, {MPSolver::CPLEX_LINEAR_PROGRAMMING, "cplex_lp"}, + {MPSolver::MOSEK_LINEAR_PROGRAMMING, "mosek_lp"}, {MPSolver::XPRESS_LINEAR_PROGRAMMING, "xpress_lp"}, {MPSolver::SCIP_MIXED_INTEGER_PROGRAMMING, "scip"}, {MPSolver::CBC_MIXED_INTEGER_PROGRAMMING, "cbc"}, @@ -584,6 +602,7 @@ constexpr {MPSolver::PDLP_LINEAR_PROGRAMMING, "pdlp"}, {MPSolver::KNAPSACK_MIXED_INTEGER_PROGRAMMING, "knapsack"}, {MPSolver::CPLEX_MIXED_INTEGER_PROGRAMMING, "cplex"}, + {MPSolver::MOSEK_MIXED_INTEGER_PROGRAMMING, "mosek"}, {MPSolver::XPRESS_MIXED_INTEGER_PROGRAMMING, "xpress"}, }; // static diff --git a/ortools/linear_solver/linear_solver.h b/ortools/linear_solver/linear_solver.h index 0b363d5d307..65256385ba0 100644 --- a/ortools/linear_solver/linear_solver.h +++ b/ortools/linear_solver/linear_solver.h @@ -237,6 +237,8 @@ class MPSolver { XPRESS_MIXED_INTEGER_PROGRAMMING = 102, COPT_LINEAR_PROGRAMMING = 103, COPT_MIXED_INTEGER_PROGRAMMING = 104, + MOSEK_LINEAR_PROGRAMMING = 203, + MOSEK_MIXED_INTEGER_PROGRAMMING = 204, }; /// Create a solver with the given name and underlying solver backend. @@ -273,6 +275,8 @@ class MPSolver { * - GUROBI_MIXED_INTEGER_PROGRAMMING or GUROBI or GUROBI_MIP * - CPLEX_LINEAR_PROGRAMMING or CPLEX_LP * - CPLEX_MIXED_INTEGER_PROGRAMMING or CPLEX or CPLEX_MIP + * - MOSEK_LINEAR_PROGRAMMING or MOSEK_LP + * - MOSEK_MIXED_INTEGER_PROGRAMMING or MOSEK or MOSEK_MIP * - XPRESS_LINEAR_PROGRAMMING or XPRESS_LP * - XPRESS_MIXED_INTEGER_PROGRAMMING or XPRESS or XPRESS_MIP * - GLPK_LINEAR_PROGRAMMING or GLPK_LP @@ -889,6 +893,7 @@ class MPSolver { friend class CBCInterface; friend class SCIPInterface; friend class GurobiInterface; + friend class MosekInterface; friend class CplexInterface; friend class XpressInterface; friend class SLMInterface; @@ -1236,6 +1241,7 @@ class MPVariable { friend class SLMInterface; friend class GurobiInterface; friend class CplexInterface; + friend class MosekInterface; friend class XpressInterface; friend class GLOPInterface; friend class MPVariableSolutionValueTest; @@ -1382,6 +1388,7 @@ class MPConstraint { friend class CLPInterface; friend class GLPKInterface; friend class SCIPInterface; + friend class MosekInterface; friend class SLMInterface; friend class GurobiInterface; friend class CplexInterface; diff --git a/ortools/linear_solver/linear_solver.proto b/ortools/linear_solver/linear_solver.proto index d90e4a074ac..dcdded0ba57 100644 --- a/ortools/linear_solver/linear_solver.proto +++ b/ortools/linear_solver/linear_solver.proto @@ -424,6 +424,7 @@ message MPModelRequest { XPRESS_LINEAR_PROGRAMMING = 101; // Commercial, needs a valid license. NOLINT CPLEX_LINEAR_PROGRAMMING = 10; // Commercial, needs a valid license. NOLINT + MOSEK_LINEAR_PROGRAMMING = 110; // Commercial, needs a valid license. NOLINT HIGHS_LINEAR_PROGRAMMING = 15; SCIP_MIXED_INTEGER_PROGRAMMING = 3; // Recommended default for MIP models. @@ -434,6 +435,8 @@ message MPModelRequest { 102; // Commercial, needs a valid license. NOLINT CPLEX_MIXED_INTEGER_PROGRAMMING = 11; // Commercial, needs a valid license. NOLINT + MOSEK_MIXED_INTEGER_PROGRAMMING = + 112; // Commercial, needs a valid license. NOLINT HIGHS_MIXED_INTEGER_PROGRAMMING = 16; BOP_INTEGER_PROGRAMMING = 12; diff --git a/ortools/linear_solver/mosek_interface.cc b/ortools/linear_solver/mosek_interface.cc new file mode 100644 index 00000000000..5c69c0496ff --- /dev/null +++ b/ortools/linear_solver/mosek_interface.cc @@ -0,0 +1,1402 @@ +// Copyright 2010-2024 Google LLg +// 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. + +// Mosek backend to MPSolver. +// +#include + +#include +#if defined(USE_MOSEK) + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "absl/flags/flag.h" +#include "absl/log/check.h" +#include "absl/log/die_if_null.h" +#include "absl/status/status.h" +#include "absl/strings/str_format.h" +#include "absl/cleanup/cleanup.h" +#include "absl/synchronization/mutex.h" +#include "absl/time/time.h" +#include "mosek.h" +#include "ortools/base/logging.h" +#include "ortools/base/timer.h" +#include "ortools/linear_solver/linear_solver.h" +#include "ortools/linear_solver/linear_solver_callback.h" +#include "ortools/linear_solver/proto_solver/proto_utils.h" +#include "ortools/util/lazy_mutable_copy.h" +#include "ortools/util/time_limit.h" + +namespace operations_research { + +class MosekInterface : public MPSolverInterface { + public: + explicit MosekInterface(MPSolver* solver, bool mip); + ~MosekInterface() override; + + void SetOptimizationDirection(bool maximize) override; + + // ----- Solve ----- + // Solves the problem using the parameter values specified. + MPSolver::ResultStatus Solve(const MPSolverParameters& param) override; + + // ----- Directly solve ----- + // Mosek should support being interrupted, but for now we'll only support + // non-interupted solves. + bool SupportsDirectlySolveProto(std::atomic* interrupt) const override { + return false; + // return interrupt == nullptr; + } + + // Solve model dirctly, bypassing protobuffers + // MPSolutionResponse DirectlySolveProto(LazyMutableCopy + // request, + // std::atomic* interrupt) + // override; + + // Writes the model. + void Write(const std::string& filename) override; + + // ----- Model modifications and extraction ----- + // Resets extracted model + void Reset() override; + + // Modifies bounds. + void SetVariableBounds(int var_index, double lb, double ub) override; + void SetVariableInteger(int var_index, bool integer) override; + void SetConstraintBounds(int row_index, double lb, double ub) override; + + // Adds Constraint incrementally. + void AddRowConstraint(MPConstraint* ct) override; + bool AddIndicatorConstraint(MPConstraint* ct) override; + // Adds variable incrementally. + void AddVariable(MPVariable* var) override; + // Changes a coefficient in a constraint. + void SetCoefficient(MPConstraint* constraint, const MPVariable* variable, + double new_value, double old_value) override; + // Clears a constraint from all its terms. + void ClearConstraint(MPConstraint* constraint) override; + // Changes a coefficient in the linear objective + void SetObjectiveCoefficient(const MPVariable* variable, + double coefficient) override; + // Changes the constant term in the linear objective. + void SetObjectiveOffset(double value) override; + // Clears the objective from all its terms. + void ClearObjective() override; + + // ------ Query statistics on the solution and the solve ------ + // Number of simplex or interior-point iterations + int64_t iterations() const override; + // Number of branch-and-bound nodes. Only available for discrete problems. + int64_t nodes() const override; + + // Returns the basis status of a row. + MPSolver::BasisStatus row_status(int constraint_index) const override; + // Returns the basis status of a column. + MPSolver::BasisStatus column_status(int variable_index) const override; + + // ----- Misc ----- + // Queries problem type. + bool IsLP() const override { return !mip_; } + bool IsMIP() const override { return mip_; } + bool IsContinuous() const override { return IsLP(); } + + void ExtractNewVariables() override; + void ExtractNewConstraints() override; + void ExtractObjective() override; + + int SolutionCount(); + std::string SolverVersion() const override { + int major, minor, rev; + + MSK_getversion(&major, &minor, &rev); + return absl::StrFormat("Mosek library version %d.%d.%d\n", major, minor, + rev); + } + + bool InterruptSolve() override { + break_solver_ = true; + return true; + } + + void* underlying_solver() override { return reinterpret_cast(task_); } + + double ComputeExactConditionNumber() const override { + if (!IsContinuous()) { + LOG(DFATAL) << "ComputeExactConditionNumber not implemented for" + << " MOSEK_MIXED_INTEGER_PROGRAMMING"; + return 0.0; + } + + LOG(DFATAL) << "ComputeExactConditionNumber not implemented for" + << " MOSEK_LINEAR_PROGRAMMING"; + return 0.0; + } + + // Iterates through the solutions in Mosek's solution pool. + bool NextSolution() override; + + void SetCallback(MPCallback* mp_callback) override; + bool SupportsCallbacks() const override { return true; } + + void SetParameters(const MPSolverParameters&) override; + bool SetSolverSpecificParametersAsString( + const std::string& parameters) override; + + void SetRelativeMipGap(double value) override; + void SetPrimalTolerance(double value) override; + void SetDualTolerance(double value) override; + void SetPresolveMode(int value) override; + void SetScalingMode(int value) override; + void SetLpAlgorithm(int value) override; + + private: + static MSKboundkeye bk_from_bounds(double lb, double ub); + + void CheckedMosekCall(MSKrescodee r) const; + + //----- Variables ----- + + // The underlying MOSEK task, which is kept updated with all Add* calls. + MSKtask_t task_; + bool break_solver_; + std::unique_ptr ptask_; + + bool mip_; + // int current_solution_index_; + // mp callback function + MPCallback* callback_ = nullptr; + // Has length equal to the number of MPConstraints in + // MPSolverInterface::solver_. Non-negative values are indexes of the + // corresponding linear constraint in Mosek, negative index i means + // disjunctive constraint (-i-1), used for indicator constraints. + std::vector mp_cons_to_mosek_cons_; + std::vector indcon_afeidx; + + int64_t domidx_rfree; + int64_t domidx_rzero; + int64_t domidx_rplus; + int64_t domidx_rminus; +}; // MosekInterface + +namespace { + +void MosekCloneParameters(MSKtask_t tdst, MSKtask_t tsrc) { + for (int p = MSK_DPAR_BEGIN; p < MSK_DPAR_END; ++p) { + double parval; + MSK_getdouparam(tsrc, (MSKdparame)p, &parval); + MSK_putdouparam(tdst, (MSKdparame)p, parval); + } + for (int p = MSK_IPAR_BEGIN; p < MSK_IPAR_END; ++p) { + int parval; + MSK_getintparam(tsrc, (MSKiparame)p, &parval); + MSK_putintparam(tdst, (MSKiparame)p, parval); + } +} + +std::pair MosekLastError(MSKtask_t task) { + int64_t lastmsgsize; + std::vector lastmsg; + MSKrescodee lastr; + MSKrescodee r = MSK_getlasterror64(task, &lastr, 0, &lastmsgsize, nullptr); + if (MSK_RES_OK == r) { + lastmsg.resize(lastmsgsize + 1); + MSK_getlasterror64(task, &lastr, lastmsgsize, &lastmsgsize, lastmsg.data()); + return std::make_pair(lastmsg.data(), lastr); + } + + return std::make_pair("", MSK_RES_OK); +} + +void CheckedMosekCall(MSKtask_t task, MSKrescodee r) { + CHECK_EQ(MSK_RES_OK, r) << "Mosek Error " << r << ": " + << MosekLastError(task).second; +} + +// This class provides a means of interacting with the task from the callback +// function. +class MosekMPCallbackContext : public MPCallbackContext { + public: + MosekMPCallbackContext(MSKtask_t task); + + // Implementation of the interface. + MPCallbackEvent Event() override; + bool CanQueryVariableValues() override; + double VariableValue(const MPVariable* variable) override; + void AddCut(const LinearRange& cutting_plane) override; + void AddLazyConstraint(const LinearRange& lazy_constraint) override; + double SuggestSolution( + const absl::flat_hash_map& solution) override; + int64_t NumExploredNodes() override; + + void Update(MSKcallbackcodee where, const double* dinf_, const int* iinf_, + const int64_t* liinf_); + void Update(const char* msg_) { + msg = msg_; + ev = MPCallbackEvent::kMessage; + } + void Reset(); + + private: + MSKtask_t task; + // current event + MPCallbackEvent ev; + // current message if the current event is kMessage + const char* msg; + std::vector mosek_variable_values_; + // NOTE: information items are assigned in callbacks and are valid for the + // diration of that callback. + const double* dinf; + const int* iinf; + const int64_t* liinf; +}; + +void MosekMPCallbackContext::Reset() { + dinf = nullptr; + iinf = nullptr; + liinf = nullptr; + msg = nullptr; +} +void MosekMPCallbackContext::Update(MSKcallbackcodee where, const double* dinf_, + const int* iinf_, const int64_t* liinf_) { + dinf = dinf_; + iinf = iinf_; + liinf = liinf_; + // kUnknown, + // // For regaining control of the main thread in single threaded + // applications, + // // not for interacting with the solver. + // kPolling, + // // The solver is currently running presolve. + // kPresolve, + // // The solver is currently running the simplex method. + // kSimplex, + // // The solver is in the MIP loop (called periodically before starting a + // new + // // node). Useful to early termination. + // kMip, + // // Called every time a new MIP incumbent is found. + // kMipSolution, + // // Called once per pass of the cut loop inside each MIP node. + // kMipNode, + // // Called in each iterate of IPM/barrier method. + // kBarrier, + // // The solver is about to log out a message, use this callback to capture + // it. kMessage, + // // The solver is in multi-objective optimization. + // kMultiObj, + + switch (where) { + // The callback function is called after a new integer solution has been + // located by the mixed-integer optimizer. + case MSK_CALLBACK_NEW_INT_MIO: + MSK_getxx(task, MSK_SOL_ITG, mosek_variable_values_.data()); + ev = MPCallbackEvent::kMipSolution; + break; + + case MSK_CALLBACK_BEGIN_DUAL_SIMPLEX: + case MSK_CALLBACK_BEGIN_DUAL_SIMPLEX_BI: + case MSK_CALLBACK_BEGIN_PRIMAL_SIMPLEX: + case MSK_CALLBACK_BEGIN_PRIMAL_SIMPLEX_BI: + case MSK_CALLBACK_END_DUAL_SIMPLEX: + case MSK_CALLBACK_END_DUAL_SIMPLEX_BI: + case MSK_CALLBACK_END_PRIMAL_SIMPLEX: + case MSK_CALLBACK_END_PRIMAL_SIMPLEX_BI: + case MSK_CALLBACK_IM_PRIMAL_SIMPLEX: + case MSK_CALLBACK_IM_SIMPLEX: + case MSK_CALLBACK_IM_SIMPLEX_BI: + case MSK_CALLBACK_PRIMAL_SIMPLEX: + case MSK_CALLBACK_UPDATE_DUAL_SIMPLEX: + case MSK_CALLBACK_UPDATE_DUAL_SIMPLEX_BI: + case MSK_CALLBACK_UPDATE_PRIMAL_SIMPLEX: + case MSK_CALLBACK_UPDATE_PRIMAL_SIMPLEX_BI: + case MSK_CALLBACK_UPDATE_SIMPLEX: + ev = MPCallbackEvent::kSimplex; + break; + + case MSK_CALLBACK_INTPNT: + case MSK_CALLBACK_BEGIN_CONIC: + ev = MPCallbackEvent::kBarrier; + break; + + case MSK_CALLBACK_BEGIN_PRESOLVE: + ev = MPCallbackEvent::kPresolve; + break; + + case MSK_CALLBACK_BEGIN_MIO: + ev = MPCallbackEvent::kMip; + break; + + // The callback function is called from within the basis identification + // procedure when the primal phase is started. + case MSK_CALLBACK_BEGIN_PRIMAL_BI: + // Begin primal feasibility repair. + case MSK_CALLBACK_BEGIN_PRIMAL_REPAIR: + // Primal sensitivity analysis is started. + case MSK_CALLBACK_BEGIN_PRIMAL_SENSITIVITY: + // The callback function is called when the primal BI setup is started. + case MSK_CALLBACK_BEGIN_PRIMAL_SETUP_BI: + // The basis identification procedure has been started. + case MSK_CALLBACK_BEGIN_BI: + // The callback function is called from within the basis identification + // procedure when the dual phase is started. + case MSK_CALLBACK_BEGIN_DUAL_BI: + // Dual sensitivity analysis is started. + case MSK_CALLBACK_BEGIN_DUAL_SENSITIVITY: + // The callback function is called when the dual BI phase is started. + case MSK_CALLBACK_BEGIN_DUAL_SETUP_BI: + // The callback function is called when the infeasibility analyzer is + // started. + case MSK_CALLBACK_BEGIN_INFEAS_ANA: + // The callback function is called when the interior-point optimizer is + // started. + case MSK_CALLBACK_BEGIN_INTPNT: + // Begin waiting for license. + case MSK_CALLBACK_BEGIN_LICENSE_WAIT: + // The callback function is called when the optimizer is started. + case MSK_CALLBACK_BEGIN_OPTIMIZER: + // Begin QCQO reformulation. + case MSK_CALLBACK_BEGIN_QCQO_REFORMULATE: + // The callback function is called when root cut generation is started. + case MSK_CALLBACK_BEGIN_ROOT_CUTGEN: + // The callback function is called when the simplex optimizer is started. + case MSK_CALLBACK_BEGIN_SIMPLEX: + // The callback function is called from within the basis identification + // procedure when the simplex clean-up phase is started. + case MSK_CALLBACK_BEGIN_SIMPLEX_BI: + // The callback function is called when solution of root relaxation is + // started. + case MSK_CALLBACK_BEGIN_SOLVE_ROOT_RELAX: + // Begin conic reformulation. + case MSK_CALLBACK_BEGIN_TO_CONIC: + // The callback function is called from within the conic optimizer after the + // information database has been updated. + case MSK_CALLBACK_CONIC: + // The callback function is called from within the dual simplex optimizer. + case MSK_CALLBACK_DUAL_SIMPLEX: + // The callback function is called when the basis identification procedure + // is terminated. + case MSK_CALLBACK_END_BI: + // The callback function is called when the conic optimizer is terminated. + case MSK_CALLBACK_END_CONIC: + // The callback function is called from within the basis identification + // procedure when the dual phase is terminated. + case MSK_CALLBACK_END_DUAL_BI: + // Dual sensitivity analysis is terminated. + case MSK_CALLBACK_END_DUAL_SENSITIVITY: + // The callback function is called when the dual BI phase is terminated. + case MSK_CALLBACK_END_DUAL_SETUP_BI: + // The callback function is called when the infeasibility analyzer is + // terminated. + case MSK_CALLBACK_END_INFEAS_ANA: + // The callback function is called when the interior-point optimizer is + // terminated. + case MSK_CALLBACK_END_INTPNT: + // End waiting for license. + case MSK_CALLBACK_END_LICENSE_WAIT: + // The callback function is called when the mixed-integer optimizer is + // terminated. + case MSK_CALLBACK_END_MIO: + // The callback function is called when the optimizer is terminated. + case MSK_CALLBACK_END_OPTIMIZER: + // The callback function is called when the presolve is completed. + case MSK_CALLBACK_END_PRESOLVE: + // The callback function is called from within the basis identification + // procedure when the primal phase is terminated. + case MSK_CALLBACK_END_PRIMAL_BI: + // End primal feasibility repair. + case MSK_CALLBACK_END_PRIMAL_REPAIR: + // Primal sensitivity analysis is terminated. + case MSK_CALLBACK_END_PRIMAL_SENSITIVITY: + // The callback function is called when the primal BI setup is terminated. + case MSK_CALLBACK_END_PRIMAL_SETUP_BI: + // End QCQO reformulation. + case MSK_CALLBACK_END_QCQO_REFORMULATE: + // The callback function is called when root cut generation is terminated. + case MSK_CALLBACK_END_ROOT_CUTGEN: + // The callback function is called when the simplex optimizer is terminated. + case MSK_CALLBACK_END_SIMPLEX: + // The callback function is called from within the basis identification + // procedure when the simplex clean-up phase is terminated. + case MSK_CALLBACK_END_SIMPLEX_BI: + // The callback function is called when solution of root relaxation is + // terminated. + case MSK_CALLBACK_END_SOLVE_ROOT_RELAX: + // End conic reformulation. + case MSK_CALLBACK_END_TO_CONIC: + // The callback function is called from within the basis identification + // procedure at an intermediate point. + case MSK_CALLBACK_IM_BI: + // The callback function is called at an intermediate stage within the conic + // optimizer where the information database has not been updated. + case MSK_CALLBACK_IM_CONIC: + // The callback function is called from within the basis identification + // procedure at an intermediate point in the dual phase. + case MSK_CALLBACK_IM_DUAL_BI: + // The callback function is called at an intermediate stage of the dual + // sensitivity analysis. + case MSK_CALLBACK_IM_DUAL_SENSIVITY: + // The callback function is called at an intermediate point in the dual + // simplex optimizer. + case MSK_CALLBACK_IM_DUAL_SIMPLEX: + // The callback function is called at an intermediate stage within the + // interior-point optimizer where the information database has not been + // updated. + case MSK_CALLBACK_IM_INTPNT: + // MOSEK is waiting for a license. + case MSK_CALLBACK_IM_LICENSE_WAIT: + // The callback function is called from within the LU factorization + // procedure at an intermediate point. + case MSK_CALLBACK_IM_LU: + // The callback function is called at an intermediate point in the + // mixed-integer optimizer. + case MSK_CALLBACK_IM_MIO: + // The callback function is called at an intermediate point in the + // mixed-integer optimizer while running the dual simplex optimizer. + case MSK_CALLBACK_IM_MIO_DUAL_SIMPLEX: + // The callback function is called at an intermediate point in the + // mixed-integer optimizer while running the interior-point optimizer. + case MSK_CALLBACK_IM_MIO_INTPNT: + // The callback function is called at an intermediate point in the + // mixed-integer optimizer while running the primal simplex optimizer. + case MSK_CALLBACK_IM_MIO_PRIMAL_SIMPLEX: + // The callback function is called from within the matrix ordering procedure + // at an intermediate point. + case MSK_CALLBACK_IM_ORDER: + // The callback function is called from within the presolve procedure at an + // intermediate stage. + case MSK_CALLBACK_IM_PRESOLVE: + // The callback function is called from within the basis identification + // procedure at an intermediate point in the primal phase. + case MSK_CALLBACK_IM_PRIMAL_BI: + // The callback function is called at an intermediate stage of the primal + // sensitivity analysis. + case MSK_CALLBACK_IM_PRIMAL_SENSIVITY: + // The callback function is called at an intermediate stage of the conic + // quadratic reformulation. + case MSK_CALLBACK_IM_QO_REFORMULATE: + // The callback is called from within root cut generation at an intermediate + // stage. + case MSK_CALLBACK_IM_ROOT_CUTGEN: + // The callback function is called when the mixed-integer optimizer is + // restarted. + case MSK_CALLBACK_RESTART_MIO: + // The callback function is called while the task is being solved on a + // remote server. + case MSK_CALLBACK_SOLVING_REMOTE: + // The callback function is called from within the basis identification + // procedure at an intermediate point in the dual phase. + case MSK_CALLBACK_UPDATE_DUAL_BI: + // The callback function is called from within the presolve procedure. + case MSK_CALLBACK_UPDATE_PRESOLVE: + // The callback function is called from within the basis identification + // procedure at an intermediate point in the primal phase. + case MSK_CALLBACK_UPDATE_PRIMAL_BI: + ev = MPCallbackEvent::kPolling; + break; + case MSK_CALLBACK_BEGIN_READ: + case MSK_CALLBACK_BEGIN_WRITE: + case MSK_CALLBACK_END_READ: + case MSK_CALLBACK_END_WRITE: + case MSK_CALLBACK_IM_READ: + case MSK_CALLBACK_READ_OPF: + case MSK_CALLBACK_READ_OPF_SECTION: + case MSK_CALLBACK_WRITE_OPF: + ev = MPCallbackEvent::kUnknown; + break; + } +} + +MosekMPCallbackContext::MosekMPCallbackContext(MSKtask_t task) + : task(task), ev{} { + int numvar; + MSK_getnumvar(task, &numvar); + mosek_variable_values_.resize(numvar); +} + +int64_t MosekMPCallbackContext::NumExploredNodes() { + int nnodes; + MSK_getintinf(task, MSK_IINF_MIO_NUM_SOLVED_NODES, &nnodes); + + return nnodes; +} + +MPCallbackEvent MosekMPCallbackContext::Event() { return ev; } + +bool MosekMPCallbackContext::CanQueryVariableValues() { + const MPCallbackEvent where = Event(); + if (where == MPCallbackEvent::kMipSolution) { + return true; + } + return false; +} + +double MosekMPCallbackContext::VariableValue(const MPVariable* variable) { + CHECK(variable != nullptr); + CHECK(ev == MPCallbackEvent::kMipSolution || ev == MPCallbackEvent::kMipNode) + << "You can only call VariableValue at " + << ToString(MPCallbackEvent::kMipSolution) << " or " + << ToString(MPCallbackEvent::kMipNode) + << " but called from: " << ToString(ev); + return mosek_variable_values_[variable->index()]; +} + +void MosekMPCallbackContext::AddCut(const LinearRange& cutting_plane) {} +void MosekMPCallbackContext::AddLazyConstraint( + const LinearRange& lazy_constraint) {} + +double MosekMPCallbackContext::SuggestSolution( + const absl::flat_hash_map& solution) { + return 0; +} + +struct MPCallbackWithMosekContext { + MosekMPCallbackContext* context; + MPCallback* callback; + bool* break_solver; +}; + +void MSKAPI StreamCallbackImpl(MSKuserhandle_t h, const char* msg) { + //std::cerr << msg << std::flush; + + MPCallbackWithMosekContext* const callback_with_context = + static_cast(h); + + CHECK(callback_with_context != nullptr); + CHECK(callback_with_context->context != nullptr); + callback_with_context->context->Update(msg); + + if (callback_with_context->callback) { + callback_with_context->callback->RunCallback( + callback_with_context->context); + } + + callback_with_context->context->Reset(); +} + +// NOTE(user): This function must have this exact API, because we are passing +// it to Mosek as a callback. +int MSKAPI CallbackImpl(MSKtask_t task, MSKuserhandle_t h, + MSKcallbackcodee where, const double* dinf, + const int* iinf, const int64_t* liinf) { + MPCallbackWithMosekContext* const callback_with_context = + static_cast(h); + CHECK(callback_with_context != nullptr); + CHECK(callback_with_context->context != nullptr); + + if (callback_with_context->callback) { + callback_with_context->context->Update(where, dinf, iinf, liinf); + callback_with_context->callback->RunCallback( + callback_with_context->context); + } + callback_with_context->context->Reset(); + return callback_with_context->break_solver[0] ? 1 : 0; +} + +} // namespace + +#if 0 // TODO +MPSolutionResponse MosekInterface::DirectlySolveProto(LazyMutableCopy request, + std::atomic* interrupt) override { + DCHECK_EQ(interrupt, nullptr); + const bool log_error = request->enable_internal_solver_output(); + + // Here we reuse the Mosek environment to support single-use license that + // forbids creating a second environment if one already exists. + return ConvertStatusOrMPSolutionResponse( + log_error, MosekSolveProto(std::move(request), global_env_)); +} +#endif + +void MosekInterface::CheckedMosekCall(MSKrescodee r) const { + ::operations_research::CheckedMosekCall(task_, r); +} + +// Creates a LP/MIP instance with the specified name and minimization objective. +MosekInterface::MosekInterface(MPSolver* const solver, bool mip) + : MPSolverInterface(solver), + task_(nullptr), + mip_(false), + ptask_(nullptr, MSK_deletetask), + break_solver_(false) { + CheckedMosekCall(MSK_makeemptytask(nullptr, &task_)); + ptask_.reset(&task_); + CheckedMosekCall(MSK_puttaskname(task_, solver_->Name().c_str())); + CheckedMosekCall(MSK_putobjsense(task_, maximize_ + ? MSK_OBJECTIVE_SENSE_MAXIMIZE + : MSK_OBJECTIVE_SENSE_MINIMIZE)); + + CheckedMosekCall(MSK_appendrzerodomain(task_, 1, &domidx_rzero)); + CheckedMosekCall(MSK_appendrplusdomain(task_, 1, &domidx_rplus)); + CheckedMosekCall(MSK_appendrminusdomain(task_, 1, &domidx_rminus)); +} + +MosekInterface::~MosekInterface() {} + +MSKboundkeye MosekInterface::bk_from_bounds(double lb, double ub) { + return (lb <= ub + ? (std::isfinite(lb) + ? (std::isfinite(ub) ? (lb < ub ? MSK_BK_RA : MSK_BK_FX) + : MSK_BK_LO) + : (std::isfinite(ub) ? MSK_BK_UP : MSK_BK_FR)) + : MSK_BK_RA); +} + +// ------ Model modifications and extraction ----- + +void MosekInterface::Reset() { + decltype(ptask_) old_taskp(std::move(ptask_)); + MSKtask_t old_task = task_; + CheckedMosekCall(MSK_makeemptytask(nullptr, &task_)); + ptask_.reset(&task_); + + mp_cons_to_mosek_cons_.clear(); + + MosekCloneParameters(task_, old_task); + + ResetExtractionInformation(); +} + +void MosekInterface::SetOptimizationDirection(bool maximize) { + InvalidateSolutionSynchronization(); + CheckedMosekCall(MSK_putobjsense(task_, maximize + ? MSK_OBJECTIVE_SENSE_MAXIMIZE + : MSK_OBJECTIVE_SENSE_MINIMIZE)); +} + +void MosekInterface::SetVariableBounds(int var_index, double lb, double ub) { + InvalidateSolutionSynchronization(); + MSKboundkeye bk = bk_from_bounds(lb, ub); + CheckedMosekCall(MSK_putvarbound(task_, var_index, bk, lb, ub)); +} + +void MosekInterface::SetVariableInteger(int index, bool integer) { + InvalidateSolutionSynchronization(); + + CheckedMosekCall(MSK_putvartype( + task_, index, integer ? MSK_VAR_TYPE_INT : MSK_VAR_TYPE_CONT)); +} + +void MosekInterface::SetConstraintBounds(int index, double lb, double ub) { + InvalidateSolutionSynchronization(); + if (mp_cons_to_mosek_cons_[index] >= 0) { + MSKboundkeye bk = bk_from_bounds(lb, ub); + CheckedMosekCall(MSK_putvarbound(task_, index, bk, lb, ub)); + } else { + int64_t djci = -mp_cons_to_mosek_cons_[index] - 1; + int64_t afei = indcon_afeidx[djci]; + int64_t afeidxs[4] = {afei, afei, afei + 1, afei + 1}; + double b[4] = {0.0, 1.0, lb, ub}; + int64_t termsize[2]{1, 3}; + int64_t domidxs[4]{domidx_rzero, domidx_rzero, domidx_rplus, domidx_rminus}; + + if (lb <= ub && lb >= ub) { + domidxs[2] = domidx_rzero; + domidxs[3] = domidx_rfree; + } else { + if (lb < 0 && !std::isfinite(lb)) domidxs[2] = domidx_rfree; + if (ub > 0 && !std::isfinite(ub)) domidxs[3] = domidx_rfree; + } + + CheckedMosekCall( + MSK_putdjc(task_, djci, 4, domidxs, 4, afeidxs, b, 2, termsize)); + } +} + +// Ordinary linear constraint are added as ranged constraints. Indicator +// constraints are added as a disjunctive constraints with constraint lb <= Ax +// <= ub here K is a value, a range or a half-open range, and X is a binary +// variable as +// (X < 0.5) OR (lb < Ax AND Ax < ub) +// +void MosekInterface::AddRowConstraint(MPConstraint* const ct) { + int conidx; + CheckedMosekCall(MSK_getnumcon(task_, &conidx)); + CheckedMosekCall(MSK_appendcons(task_, 1)); + mp_cons_to_mosek_cons_.push_back(conidx); + + double lb = ct->lb(); + double ub = ct->ub(); + const std::string& name = ct->name(); + + MSKboundkeye bk = bk_from_bounds(lb, ub); + CheckedMosekCall(MSK_putconbound(task_, conidx, bk, lb, ub)); + std::vector cof; + cof.reserve(ct->coefficients_.size()); + std::vector subj; + subj.reserve(cof.size()); + for (const auto& it : ct->terms()) { + subj.push_back(it.first->index()); + cof.push_back(it.second); + } + if (name.size() > 0) + CheckedMosekCall(MSK_putconname(task_, conidx, name.c_str())); + + CheckedMosekCall( + MSK_putarow(task_, conidx, subj.size(), subj.data(), cof.data())); +} + +bool MosekInterface::AddIndicatorConstraint(MPConstraint* const ct) { + int64_t djci, afei; + CheckedMosekCall(MSK_getnumdjc(task_, &djci)); + CheckedMosekCall(MSK_appenddjcs(task_, 1)); + CheckedMosekCall(MSK_getnumafe(task_, &afei)); + CheckedMosekCall(MSK_appendafes(task_, 2)); + mp_cons_to_mosek_cons_.push_back(-1 - djci); + indcon_afeidx.push_back(afei); + + auto indvar = ct->indicator_variable(); + + CHECK(indvar != nullptr); + CheckedMosekCall(MSK_putvartype(task_, indvar->index(), MSK_VAR_TYPE_INT)); + + CheckedMosekCall( + MSK_putvarbound(task_, indvar->index(), MSK_BK_RA, 0.0, 1.0)); + + const std::string& name = ct->name(); + if (name.size() > 0) + CheckedMosekCall(MSK_putdjcname(task_, djci, name.c_str())); + + { + double lb = ct->lb(); + double ub = ct->ub(); + + int64_t domidxs[4] = {domidx_rzero, domidx_rzero, domidx_rplus, + domidx_rminus}; + int64_t afeidxs[4] = {afei, afei, afei + 1, afei + 1}; + double b[4] = {0.0, 1.0, lb, ub}; + int64_t termsize[2] = {1, 3}; + + if (lb <= ub && lb >= ub) { + domidxs[2] = domidx_rzero; + domidxs[3] = domidx_rfree; + } else { + if (lb < 0 && !std::isfinite(lb)) domidxs[2] = domidx_rfree; + if (ub > 0 && !std::isfinite(ub)) domidxs[3] = domidx_rfree; + } + + CheckedMosekCall( + MSK_putdjc(task_, djci, 4, domidxs, 4, afeidxs, b, 2, termsize)); + } + { + std::vector cof; + cof.reserve(ct->coefficients_.size()); + std::vector subj; + subj.reserve(cof.size()); + for (const auto& it : ct->terms()) { + subj.push_back(it.first->index()); + cof.push_back(it.second); + } + CheckedMosekCall( + MSK_putafefrow(task_, afei + 1, subj.size(), subj.data(), cof.data())); + } + { + double c = 1.0; + int indvarj = indvar->index(); + CheckedMosekCall(MSK_putafefrow(task_, afei, 1, &indvarj, &c)); + } + + return true; +} + +void MosekInterface::AddVariable(MPVariable* const var) { + int j; + CheckedMosekCall(MSK_getnumvar(task_, &j)); + CheckedMosekCall(MSK_appendvars(task_, 1)); + double lb = var->lb(); + double ub = var->ub(); + + const std::string& name = var->name(); + + MSKboundkeye bk = bk_from_bounds(lb, ub); + CheckedMosekCall(MSK_putvarbound(task_, j, bk, lb, ub)); + if (var->integer()) + CheckedMosekCall(MSK_putvartype(task_, j, MSK_VAR_TYPE_INT)); + if (name.size() > 0) CheckedMosekCall(MSK_putvarname(task_, j, name.c_str())); +} + +void MosekInterface::SetCoefficient(MPConstraint* const constraint, + const MPVariable* const variable, + double new_value, double old_value) { + InvalidateSolutionSynchronization(); + + ssize_t coni = mp_cons_to_mosek_cons_[constraint->index()]; + if (coni >= 0) { + CheckedMosekCall(MSK_putaij(task_, coni, variable->index(), new_value)); + } else { + int64_t djci = -coni - 1; + int64_t afei = indcon_afeidx[djci] + 1; + + CheckedMosekCall( + MSK_putafefentry(task_, afei, variable->index(), new_value)); + } +} + +// Question: Is an indicator constraint ever cleared? What exactly does that +// mean? +void MosekInterface::ClearConstraint(MPConstraint* const constraint) { + InvalidateSolutionSynchronization(); + // TODO: Cleanup if afe nonzeros etc? + ssize_t coni = mp_cons_to_mosek_cons_[constraint->index()]; + if (coni >= 0) { + CheckedMosekCall(MSK_putarow(task_, coni, 0, nullptr, nullptr)); + CheckedMosekCall(MSK_putconbound(task_, coni, MSK_BK_FR, 0.0, 0.0)); + } else { + int64_t djci = -coni - 1; + CheckedMosekCall( + MSK_putdjc(task_, djci, 0, nullptr, 0, nullptr, nullptr, 0, nullptr)); + } +} + +void MosekInterface::SetObjectiveCoefficient(const MPVariable* const variable, + double coefficient) { + InvalidateSolutionSynchronization(); + CheckedMosekCall(MSK_putcj(task_, variable->index(), coefficient)); +} + +void MosekInterface::SetObjectiveOffset(double value) { + InvalidateSolutionSynchronization(); + + CheckedMosekCall(MSK_putcfix(task_, value)); +} + +void MosekInterface::ClearObjective() { + InvalidateSolutionSynchronization(); + int numvar; + CheckedMosekCall(MSK_getnumvar(task_, &numvar)); + for (int i = 0; i < numvar; ++i) MSK_putcj(task_, i, 0.0); + MSK_putcfix(task_, 0.0); +} + +// ------ Query statistics on the solution and the solve ------ + +int64_t MosekInterface::iterations() const { + if (!CheckSolutionIsSynchronized()) return kUnknownNumberOfIterations; + // TODO: Iter Count + int32_t psim_iter, dsim_iter, intpnt_iter; + CheckedMosekCall( + MSK_getintinf(task_, MSK_IINF_SIM_DUAL_ITER, &psim_iter)); + CheckedMosekCall( + MSK_getintinf(task_, MSK_IINF_SIM_PRIMAL_ITER, &dsim_iter)); + CheckedMosekCall( + MSK_getintinf(task_, MSK_IINF_INTPNT_ITER, &intpnt_iter)); + + int64_t mio_intpnt_iter, mio_simplex_iter; + CheckedMosekCall(MSK_getlintinf(task_,MSK_LIINF_MIO_INTPNT_ITER, &mio_intpnt_iter)); + CheckedMosekCall(MSK_getlintinf(task_,MSK_LIINF_MIO_SIMPLEX_ITER, &mio_simplex_iter)); + + if (intpnt_iter > 0) + return intpnt_iter; + else if (psim_iter+dsim_iter > 0) + return psim_iter+dsim_iter; + else if (mio_simplex_iter > 0) + return mio_simplex_iter; + else if (mio_intpnt_iter > 0) + return mio_intpnt_iter; + else + return 0; +} + +int64_t MosekInterface::nodes() const { + if (!CheckSolutionIsSynchronized()) return kUnknownNumberOfNodes; + int nnodes; + CheckedMosekCall( + MSK_getintinf(task_, MSK_IINF_MIO_NUM_SOLVED_NODES, &nnodes)); + return nnodes; +} + +// Returns the basis status of a row. +MPSolver::BasisStatus MosekInterface::row_status(int constraint_index) const { + auto coni = mp_cons_to_mosek_cons_[constraint_index]; + if (coni < 0) { + LOG(DFATAL) << "Basis status only available for continuous problems."; + } + + int soldef; + + CheckedMosekCall(MSK_solutiondef(task_, MSK_SOL_BAS, &soldef)); + if (!soldef) { + LOG(DFATAL) + << "Basis status only available when a basis solution has been found."; + return MPSolver::FREE; + } + + MSKstakeye sk; + CheckedMosekCall(MSK_getskcslice(task_, MSK_SOL_BAS, coni, coni + 1, &sk)); + + switch (sk) { + case MSK_SK_BAS: + return MPSolver::BASIC; + case MSK_SK_LOW: + return MPSolver::AT_LOWER_BOUND; + case MSK_SK_UPR: + return MPSolver::AT_UPPER_BOUND; + case MSK_SK_FIX: + return MPSolver::FIXED_VALUE; + case MSK_SK_SUPBAS: + return MPSolver::FREE; + + default: + LOG(DFATAL) << "Basis status only available when a basis solution has " + "been found."; + return MPSolver::FREE; + } +} + +// Returns the basis status of a column. +MPSolver::BasisStatus MosekInterface::column_status(int variable_index) const { + int soldef; + MSKsoltypee whichsol; + + if (MSK_RES_OK == MSK_solutiondef(task_, MSK_SOL_ITG, &soldef) && soldef) { + whichsol = MSK_SOL_ITG; + } else if (MSK_RES_OK == MSK_solutiondef(task_, MSK_SOL_BAS, &soldef) && + soldef) { + whichsol = MSK_SOL_BAS; + } else if (MSK_RES_OK == MSK_solutiondef(task_, MSK_SOL_ITR, &soldef) && + soldef) { + whichsol = MSK_SOL_ITG; + } + + if (!soldef) { + LOG(DFATAL) + << "Basis status only available when a basis solution has been found."; + return MPSolver::FREE; + } + + MSKstakeye sk; + CheckedMosekCall(MSK_getskxslice(task_, whichsol, variable_index, + variable_index + 1, &sk)); + + switch (sk) { + case MSK_SK_BAS: + return MPSolver::BASIC; + case MSK_SK_LOW: + return MPSolver::AT_LOWER_BOUND; + case MSK_SK_UPR: + return MPSolver::AT_UPPER_BOUND; + case MSK_SK_FIX: + return MPSolver::FIXED_VALUE; + case MSK_SK_SUPBAS: + return MPSolver::FREE; + default: + return MPSolver::FREE; + } +} + +// Extracts new variables. +void MosekInterface::ExtractNewVariables() { + int numvar; + auto total_num_vars = solver_->variables_.size(); + CheckedMosekCall(MSK_getnumvar(task_, &numvar)); + if (total_num_vars > numvar) { + CheckedMosekCall(MSK_appendvars(task_, total_num_vars - numvar)); + auto& obj = solver_->Objective(); + + for (int j = numvar; numvar < total_num_vars; ++j) { + auto var = solver_->variables_[j]; + set_variable_as_extracted(j, true); + MSKboundkeye bk = bk_from_bounds(var->lb(), var->ub()); + CheckedMosekCall(MSK_putvarbound(task_, j, bk, var->lb(), var->ub())); + if (var->integer()) { + CheckedMosekCall(MSK_putvartype(task_, j, MSK_VAR_TYPE_INT)); + } + + double cj = obj.GetCoefficient(var); + if (cj > 0 || cj < 0) { + CheckedMosekCall(MSK_putcj(task_, j, cj)); + } + + for (int i = 0; i < mp_cons_to_mosek_cons_.size(); ++i) { + auto coni = mp_cons_to_mosek_cons_[i]; + const MPConstraint* ct = solver_->constraints()[i]; + if (coni >= 0) { + for (const auto& item : ct->coefficients_) { + if (item.first->index() >= numvar) { + CheckedMosekCall( + MSK_putaij(task_, coni, item.first->index(), item.second)); + } + } + } else { + auto djci = -coni - 1; + auto afei = indcon_afeidx[djci]; + for (const auto& item : ct->coefficients_) { + if (item.first->index() >= numvar) { + CheckedMosekCall(MSK_putafefentry( + task_, afei + 1, item.first->index(), item.second)); + } + } + } + } + } + } +} + +void MosekInterface::ExtractNewConstraints() { + int total_num_rows = solver_->constraints_.size(); + if (mp_cons_to_mosek_cons_.size() < total_num_rows) { + // Add each new constraint. + for (int row = last_constraint_index_; row < total_num_rows; ++row) { + MPConstraint* const ct = solver_->constraints_[row]; + set_constraint_as_extracted(row, true); + AddRowConstraint(ct); + } + } +} + +void MosekInterface::ExtractObjective() { + CheckedMosekCall(MSK_putobjsense(task_, maximize_ + ? MSK_OBJECTIVE_SENSE_MAXIMIZE + : MSK_OBJECTIVE_SENSE_MINIMIZE)); + const auto& obj = solver_->Objective(); + CheckedMosekCall(MSK_putcfix(task_, obj.offset())); +} + +// ------ Parameters ----- + +void MosekInterface::SetParameters(const MPSolverParameters& param) { + SetCommonParameters(param); + if (mip_) { + SetMIPParameters(param); + } +} + +bool MosekInterface::SetSolverSpecificParametersAsString( + const std::string& parameters) { + std::string_view data(parameters); + std::string key, value; + while (data.size() > 0) { + auto p = data.find('\n'); + if (p == std::string_view::npos) p = data.size(); + std::string_view line = data.substr(0, p); + + auto eq_pos = line.find('='); + if (eq_pos != std::string_view::npos) { + key.clear(); + value.clear(); + key += line.substr(0, eq_pos); + value += line.substr(eq_pos + 1); + + if (MSK_RES_OK != MSK_putparam(task_, key.c_str(), value.c_str())) { + LOG(WARNING) << "Failed to set parameters '" << key << "' to '" << value + << "'"; + } + } + + data = data.substr(p + 1); + } + return true; +} + +void MosekInterface::SetRelativeMipGap(double value) { + CheckedMosekCall(MSK_putdouparam(task_, MSK_DPAR_MIO_REL_GAP_CONST, value)); +} + +void MosekInterface::SetPrimalTolerance(double value) { + CheckedMosekCall(MSK_putdouparam(task_, MSK_DPAR_INTPNT_TOL_PFEAS, value)); + CheckedMosekCall(MSK_putdouparam(task_, MSK_DPAR_BASIS_TOL_X, value)); +} + +void MosekInterface::SetDualTolerance(double value) { + CheckedMosekCall(MSK_putdouparam(task_, MSK_DPAR_INTPNT_TOL_DFEAS, value)); + CheckedMosekCall(MSK_putdouparam(task_, MSK_DPAR_BASIS_TOL_S, value)); +} + +void MosekInterface::SetPresolveMode(int value) { + switch (value) { + case MPSolverParameters::PRESOLVE_OFF: { + CheckedMosekCall(MSK_putintparam(task_, MSK_IPAR_PRESOLVE_USE, MSK_OFF)); + break; + } + case MPSolverParameters::PRESOLVE_ON: { + CheckedMosekCall(MSK_putintparam(task_, MSK_IPAR_PRESOLVE_USE, MSK_ON)); + break; + } + default: { + SetIntegerParamToUnsupportedValue(MPSolverParameters::PRESOLVE, value); + } + } +} + +// Sets the scaling mode. +void MosekInterface::SetScalingMode(int value) { + switch (value) { + case MPSolverParameters::SCALING_OFF: + CheckedMosekCall( + MSK_putintparam(task_, MSK_IPAR_INTPNT_SCALING, MSK_SCALING_NONE)); + CheckedMosekCall( + MSK_putintparam(task_, MSK_IPAR_SIM_SCALING, MSK_SCALING_NONE)); + break; + case MPSolverParameters::SCALING_ON: + CheckedMosekCall( + MSK_putintparam(task_, MSK_IPAR_INTPNT_SCALING, MSK_SCALING_FREE)); + CheckedMosekCall( + MSK_putintparam(task_, MSK_IPAR_SIM_SCALING, MSK_SCALING_FREE)); + break; + default: + // Leave the parameters untouched. + break; + } +} + +void MosekInterface::SetLpAlgorithm(int value) { + switch (value) { + case MPSolverParameters::DUAL: + CheckedMosekCall(MSK_putintparam(task_, MSK_IPAR_OPTIMIZER, + MSK_OPTIMIZER_DUAL_SIMPLEX)); + break; + case MPSolverParameters::PRIMAL: + CheckedMosekCall(MSK_putintparam(task_, MSK_IPAR_OPTIMIZER, + MSK_OPTIMIZER_PRIMAL_SIMPLEX)); + break; + case MPSolverParameters::BARRIER: + CheckedMosekCall( + MSK_putintparam(task_, MSK_IPAR_OPTIMIZER, MSK_OPTIMIZER_INTPNT)); + break; + default: + SetIntegerParamToUnsupportedValue(MPSolverParameters::LP_ALGORITHM, + value); + } +} + +int MosekInterface::SolutionCount() { + int soldef; + MSK_solutiondef(task_, MSK_SOL_ITG, &soldef); + if (soldef) return 1; + MSK_solutiondef(task_, MSK_SOL_BAS, &soldef); + if (soldef) return 1; + MSK_solutiondef(task_, MSK_SOL_ITR, &soldef); + if (soldef) return 1; + return 0; +} + +MPSolver::ResultStatus MosekInterface::Solve(const MPSolverParameters& param) { + WallTimer timer; + timer.Start(); + + + // Set log level. + CheckedMosekCall(MSK_putintparam(task_, MSK_IPAR_LOG, quiet_ ? 0 : 10)); + + ExtractModel(); + // Sync solver. + VLOG(1) << absl::StrFormat("Model built in %s.", + absl::FormatDuration(timer.GetDuration())); + + int numvar; + MSK_getnumvar(task_, &numvar); + bool ismip = false; + for (int j = 0; j < numvar && ! ismip; ++j) { + MSKvariabletypee vt; + MSK_getvartype(task_,j,&vt); ismip = ismip && vt == MSK_VAR_TYPE_INT; + } + // Set solution hints. Currently this only affects integer solution. + if (solver_->solution_hint_.size() > 0) { + std::vector xx(numvar); + for (const std::pair& p : + solver_->solution_hint_) { + xx[p.first->index()] = p.second; + } + MSK_putxx(task_, MSK_SOL_ITG, xx.data()); + } + + // Time limit. + if (solver_->time_limit() != 0) { + VLOG(1) << "Setting time limit = " << solver_->time_limit() << " ms."; + CheckedMosekCall(MSK_putdouparam(task_, MSK_DPAR_OPTIMIZER_MAX_TIME, + solver_->time_limit_in_secs())); + } + + // We first set our internal MPSolverParameters from 'param' and then set + // any user-specified internal solver parameters via + // solver_specific_parameter_string_. + // Default MPSolverParameters can override custom parameters (for example for + // presolving) and therefore we apply MPSolverParameters first. + SetParameters(param); + solver_->SetSolverSpecificParametersAsString( + solver_->solver_specific_parameter_string_); + + // remove any pre-existing solution in task that are not relevant for the + // result. + MSK_putintparam(task_, MSK_IPAR_REMOVE_UNUSED_SOLUTIONS, MSK_ON); + + // Solve + timer.Restart(); + + MSKrescodee trm; + { + MosekMPCallbackContext mosek_context(task_); + MPCallback* cb = callback_; + MPCallbackWithMosekContext mp_callback_with_context = { + .context = &mosek_context, + .callback = cb, + .break_solver = &break_solver_}; + + auto _remove_cb = absl::MakeCleanup([&]() { + MSK_linkfunctotaskstream(task_, MSK_STREAM_LOG, nullptr, nullptr); + MSK_putcallbackfunc(task_, nullptr, nullptr); + }); + MSK_putcallbackfunc(task_, CallbackImpl, &mp_callback_with_context); + MSK_linkfunctotaskstream(task_, MSK_STREAM_LOG, &mp_callback_with_context, + StreamCallbackImpl); + CheckedMosekCall(MSK_optimizetrm(task_, &trm)); + } + + VLOG(1) << absl::StrFormat("Solved in %s.", + absl::FormatDuration(timer.GetDuration())); + // Get the status. + MSKprostae prosta = (MSKprostae)-1; + MSKsolstae solsta = (MSKsolstae)-1; + MSKsoltypee whichsol; + bool soldefined = true; + { + int soldef; + whichsol = MSK_SOL_ITG; + MSK_solutiondef(task_, whichsol, &soldef); + if (!soldef) { + whichsol = MSK_SOL_BAS; + MSK_solutiondef(task_, whichsol, &soldef); + } + if (!soldef) { + whichsol = MSK_SOL_ITR; + MSK_solutiondef(task_, whichsol, &soldef); + } + soldefined = soldef != 0; + } + + if (soldefined) { + MSK_getprosta(task_, whichsol, &prosta); + MSK_getsolsta(task_, whichsol, &solsta); + } + + VLOG(1) << absl::StrFormat("Solution status %d.\n", prosta); + const int solution_count = SolutionCount(); + + if (!soldefined) { + result_status_ = MPSolver::NOT_SOLVED; + } else if (solsta == MSK_SOL_STA_OPTIMAL || + solsta == MSK_SOL_STA_INTEGER_OPTIMAL) { + result_status_ = MPSolver::OPTIMAL; + } else if (solsta == MSK_SOL_STA_PRIM_AND_DUAL_FEAS) { + result_status_ = MPSolver::FEASIBLE; + } else if (prosta == MSK_PRO_STA_PRIM_INFEAS) { + result_status_ = MPSolver::INFEASIBLE; + } else if (prosta == MSK_PRO_STA_DUAL_INFEAS) { + result_status_ = MPSolver::UNBOUNDED; + } else if (prosta == MSK_PRO_STA_PRIM_INFEAS_OR_UNBOUNDED) { + // TODO(user): We could introduce our own "infeasible or + // unbounded" status. + result_status_ = MPSolver::INFEASIBLE; + } else { + result_status_ = MPSolver::NOT_SOLVED; + } + + // Get best objective bound value + if (whichsol == MSK_SOL_ITG && + (result_status_ == MPSolver::FEASIBLE || + result_status_ == MPSolver::OPTIMAL)) { + MSK_getdouinf(task_, MSK_DINF_MIO_OBJ_BOUND, &best_objective_bound_); + VLOG(1) << "best bound = " << best_objective_bound_; + } + + if (solution_count > 0 && (result_status_ == MPSolver::FEASIBLE || + result_status_ == MPSolver::OPTIMAL)) { + // Get the results. + MSK_getprimalobj(task_, whichsol, &objective_value_); + VLOG(1) << "objective = " << objective_value_; + + std::vector xx(numvar); + CheckedMosekCall(MSK_getxx(task_, whichsol, xx.data())); + { + for (int i = 0; i < solver_->variables_.size(); ++i) { + MPVariable* const var = solver_->variables_[i]; + var->set_solution_value(xx[i]); + VLOG(3) << var->name() << ", value = " << xx[i]; + } + } + if (whichsol != MSK_SOL_ITG) { + { + std::vector slx(numvar); + std::vector sux(numvar); + + CheckedMosekCall(MSK_getslx(task_, whichsol, slx.data())); + CheckedMosekCall(MSK_getsux(task_, whichsol, sux.data())); + + for (int i = 0; i < solver_->variables_.size(); ++i) { + MPVariable* const var = solver_->variables_[i]; + var->set_reduced_cost(slx[i] - sux[i]); + VLOG(4) << var->name() << ", reduced cost = " << (slx[i] - sux[i]); + } + } + + { + size_t numcon = mp_cons_to_mosek_cons_.size(); + std::vector y(numcon); + + CheckedMosekCall(MSK_gety(task_, whichsol, y.data())); + + for (int i = 0; i < solver_->constraints_.size(); ++i) { + MPConstraint* const ct = solver_->constraints_[i]; + auto coni = mp_cons_to_mosek_cons_[ct->index()]; + if (coni >= 0) { + ct->set_dual_value(y[coni]); + VLOG(4) << "row " << ct->index() << ", dual value = " << y[coni]; + } + } + } + } + } + + sync_status_ = SOLUTION_SYNCHRONIZED; + return result_status_; +} + +// Select next solution and assign all solution values to variables. +bool MosekInterface::NextSolution() { return false; } + +void MosekInterface::Write(const std::string& filename) { + // if (sync_status_ == MUST_RELOAD) { + // Reset(); + // } + // ExtractModel(); + VLOG(1) << "Writing Mosek Task file \"" << filename << "\"."; + MSKrescodee r = MSK_writedata(task_, filename.c_str()); + if (MSK_RES_OK != r) { + auto lasterr = MosekLastError(task_); + LOG(WARNING) << "Failed to write Task. Error (" << lasterr.second + << "): " << lasterr.first; + } +} + +MPSolverInterface* BuildMosekInterface(bool mip, MPSolver* const solver) { + return new MosekInterface(solver, mip); +} + +void MosekInterface::SetCallback(MPCallback* mp_callback) { + callback_ = mp_callback; +} + +} // namespace operations_research + +#endif // USE_MOSEK diff --git a/ortools/math_opt/cpp/parameters.h b/ortools/math_opt/cpp/parameters.h index ba25d991226..d8d85f4765f 100644 --- a/ortools/math_opt/cpp/parameters.h +++ b/ortools/math_opt/cpp/parameters.h @@ -109,6 +109,8 @@ enum class SolverType { // Slow/not recommended for production. Not an LP solver (no dual information // returned). kSantorini = SOLVER_TYPE_SANTORINI, + + kMosek = SOLVER_TYPE_MOSEK, }; MATH_OPT_DEFINE_ENUM(SolverType, SOLVER_TYPE_UNSPECIFIED); diff --git a/ortools/math_opt/parameters.proto b/ortools/math_opt/parameters.proto index af0c01bb72c..a519f064d40 100644 --- a/ortools/math_opt/parameters.proto +++ b/ortools/math_opt/parameters.proto @@ -105,6 +105,12 @@ enum SolverTypeProto { // Slow/not recommended for production. Not an LP solver (no dual information // returned). SOLVER_TYPE_SANTORINI = 11; + + // Mosek solver (third party). + // + // Supports LP, MIP, and nonconvex integer quadratic problems. Generally the + // fastest option, but has special licensing. + SOLVER_TYPE_MOSEK = 12; } // Selects an algorithm for solving linear programs. diff --git a/ortools/math_opt/python/parameters.py b/ortools/math_opt/python/parameters.py index e3416557080..cd742ad0f32 100644 --- a/ortools/math_opt/python/parameters.py +++ b/ortools/math_opt/python/parameters.py @@ -76,6 +76,8 @@ class SolverType(enum.Enum): QPs are unimplemented). SANTORINI: The Santorini Solver (first party). Supports MIP. Experimental, do not use in production. + MOSEK: Mosek solver (third party). Supports LP, MIP, conic quadratic. + Requires a license. """ GSCIP = math_opt_parameters_pb2.SOLVER_TYPE_GSCIP @@ -89,6 +91,7 @@ class SolverType(enum.Enum): SCS = math_opt_parameters_pb2.SOLVER_TYPE_SCS HIGHS = math_opt_parameters_pb2.SOLVER_TYPE_HIGHS SANTORINI = math_opt_parameters_pb2.SOLVER_TYPE_SANTORINI + MOSEK = math_opt_parameters_pb2.SOLVER_TYPE_MOSEK def solver_type_from_proto( diff --git a/ortools/math_opt/solvers/BUILD.bazel b/ortools/math_opt/solvers/BUILD.bazel index 023e966a490..f43f21f99aa 100644 --- a/ortools/math_opt/solvers/BUILD.bazel +++ b/ortools/math_opt/solvers/BUILD.bazel @@ -102,6 +102,59 @@ cc_test( ], ) +cc_library( + name = "mosek_solver", + srcs = [ + "mosek_solver.cc", + ], + hdrs = [ + "mosek_solver.h", + ], + visibility = ["//visibility:public"], + deps = [ + "//ortools/base:linked_hash_map", + "//ortools/base:map_util", + "//ortools/base:protoutil", + "//ortools/base:status_macros", + "//ortools/gurobi:environment", + "//ortools/gurobi/isv_public:gurobi_isv", + "//ortools/math_opt:callback_cc_proto", + "//ortools/math_opt:infeasible_subsystem_cc_proto", + "//ortools/math_opt:model_cc_proto", + "//ortools/math_opt:model_parameters_cc_proto", + "//ortools/math_opt:model_update_cc_proto", + "//ortools/math_opt:parameters_cc_proto", + "//ortools/math_opt:result_cc_proto", + "//ortools/math_opt:solution_cc_proto", + "//ortools/math_opt:sparse_containers_cc_proto", + "//ortools/math_opt/core:invalid_indicators", + "//ortools/math_opt/core:inverted_bounds", + "//ortools/math_opt/core:math_opt_proto_utils", + "//ortools/math_opt/core:non_streamable_solver_init_arguments", + "//ortools/math_opt/core:solver_interface", + "//ortools/math_opt/core:sorted", + "//ortools/math_opt/core:sparse_vector_view", + "//ortools/math_opt/solvers/mosek:mosekwrp", + "//ortools/math_opt/validators:callback_validator", + "//ortools/port:proto_utils", + "//ortools/util:solve_interrupter", + "//ortools/util:testing_utils", + "@com_google_absl//absl/algorithm:container", + "@com_google_absl//absl/container:flat_hash_map", + "@com_google_absl//absl/container:flat_hash_set", + "@com_google_absl//absl/log", + "@com_google_absl//absl/log:check", + "@com_google_absl//absl/memory", + "@com_google_absl//absl/meta:type_traits", + "@com_google_absl//absl/status", + "@com_google_absl//absl/status:statusor", + "@com_google_absl//absl/strings", + "@com_google_absl//absl/time", + "@com_google_absl//absl/types:span", + "@com_google_protobuf//:protobuf", + ], +) + cc_library( name = "gurobi_callback", srcs = ["gurobi_callback.cc"], diff --git a/ortools/math_opt/solvers/CMakeLists.txt b/ortools/math_opt/solvers/CMakeLists.txt index 9cae3b1ea7f..0ab2e0ba53c 100644 --- a/ortools/math_opt/solvers/CMakeLists.txt +++ b/ortools/math_opt/solvers/CMakeLists.txt @@ -31,6 +31,11 @@ if(NOT USE_HIGHS) list(FILTER _SRCS EXCLUDE REGEX "/highs_.*.h$") list(FILTER _SRCS EXCLUDE REGEX "/highs_.*.cc$") endif() +if(NOT USE_MOSEK) + list(FILTER _SRCS EXCLUDE REGEX "/mosek/") + list(FILTER _SRCS EXCLUDE REGEX "/mosek_.*.h$") + list(FILTER _SRCS EXCLUDE REGEX "/mosek_.*.cc$") +endif() if(NOT USE_PDLP) list(FILTER _SRCS EXCLUDE REGEX "/pdlp_.*.h$") list(FILTER _SRCS EXCLUDE REGEX "/pdlp_.*.cc$") @@ -49,6 +54,7 @@ target_link_libraries(${NAME} PRIVATE absl::strings $<$:GLPK::GLPK> $<$:highs::highs> + $<$:mosek::mosek> $<$:Eigen3::Eigen> $<$:libscip> ${PROJECT_NAMESPACE}::math_opt_proto) diff --git a/ortools/math_opt/solvers/mosek.proto b/ortools/math_opt/solvers/mosek.proto new file mode 100644 index 00000000000..ed472d9953b --- /dev/null +++ b/ortools/math_opt/solvers/mosek.proto @@ -0,0 +1,28 @@ +// Copyright 2010-2024 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. + +// Proto messages specific to Mosek. +syntax = "proto3"; + +package operations_research.math_opt; + +// Mosek specific parameters for solving. See +// https://docs.mosek.com/latest/capi/parameters.html +// for a list of possible parameters. +message MosekParametersProto { + message Parameter { + string name = 1; + string value = 2; + } + repeated Parameter parameters = 1; +} diff --git a/ortools/math_opt/solvers/mosek/BUILD.bazel b/ortools/math_opt/solvers/mosek/BUILD.bazel new file mode 100644 index 00000000000..da67742fa24 --- /dev/null +++ b/ortools/math_opt/solvers/mosek/BUILD.bazel @@ -0,0 +1,39 @@ +# Copyright 2010-2024 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. + +cc_library( + name = "mosekwrp", + srcs = [ + "mosekwrp.cc", + ], + hdrs = [ + "mosekwrp.h", + ], + visibility = [ + "//ortools/math_opt:__subpackages__", + ], + deps = [ + "//ortools/base:logging", + "//ortools/base:source_location", + "//ortools/base:status_macros", + "//ortools/math_opt/solvers:mosek_cc_proto", + "@com_google_absl//absl/log:check", + "@com_google_absl//absl/log:die_if_null", + "@com_google_absl//absl/memory", + "@com_google_absl//absl/status", + "@com_google_absl//absl/status:statusor", + "@com_google_absl//absl/strings:str_format", + "@com_google_absl//absl/types:span", + ], +) + diff --git a/ortools/math_opt/solvers/mosek/mosekwrp.cc b/ortools/math_opt/solvers/mosek/mosekwrp.cc new file mode 100644 index 00000000000..f1f64420a06 --- /dev/null +++ b/ortools/math_opt/solvers/mosek/mosekwrp.cc @@ -0,0 +1,546 @@ +#include "mosekwrp.h" + +#include + +#include +#include +#include +#include +#include + +#include "absl/cleanup/cleanup.h" +#include "absl/status/status.h" + +namespace operations_research::math_opt { + +Mosek* Mosek::Create() { + MSKtask_t task; + MSKrescodee r = MSK_makeemptytask(nullptr, &task); + if (r != MSK_RES_OK) { + return nullptr; + } + int64_t domidx; + MSK_appendrzerodomain(task, 0, &domidx); + + return new Mosek(task); +} + +Mosek::Mosek(MSKtask_t task) : task(task, delete_msk_task_func) {} +Mosek::Mosek(Mosek&& m) : task(std::move(m.task)) {} + +void Mosek::PutName(const std::string& name) { + MSK_puttaskname(task.get(), name.c_str()); +} +void Mosek::PutObjName(const std::string& name) { + MSK_putobjname(task.get(), name.c_str()); +} +void Mosek::PutVarName(VariableIndex j, const std::string& name) { + MSK_putvarname(task.get(), j, name.c_str()); +} +void Mosek::PutConName(ConstraintIndex i, const std::string& name) { + MSK_putconname(task.get(), i, name.c_str()); +} +void Mosek::PutObjectiveSense(bool maximize) { + MSK_putobjsense(task.get(), maximize ? MSK_OBJECTIVE_SENSE_MAXIMIZE + : MSK_OBJECTIVE_SENSE_MINIMIZE); +} + +absl::StatusOr Mosek::AppendVars( + const std::vector& lb, const std::vector& ub) { + if (lb.size() != ub.size()) + return absl::InvalidArgumentError("Mismatching lengths of lb and ub"); + size_t n = lb.size(); + int firstj = NumVar(); + if (n > std::numeric_limits::max()) + return absl::InvalidArgumentError("arguments lb and ub too large"); + + MSK_appendvars(task.get(), (int)n); + std::vector bk(n); + for (ssize_t i = 0; i < n; ++i) { + bk[i] = (lb[i] > ub[i] + ? MSK_BK_RA + : (std::isfinite(lb[i]) + ? (std::isfinite(ub[i]) + ? (lb[i] < ub[i] ? MSK_BK_RA : MSK_BK_FX) + : MSK_BK_LO) + : (std::isfinite(ub[i]) ? MSK_BK_UP : MSK_BK_FR))); + } + + MSK_putvarboundslice(task.get(), firstj, firstj + (int)n, bk.data(), + lb.data(), ub.data()); + return firstj; +} +absl::StatusOr Mosek::AppendCons( + const std::vector& lb, const std::vector& ub) { + if (lb.size() != ub.size()) + return absl::InvalidArgumentError("Mismatching lengths of lb and ub"); + size_t n = lb.size(); + int firsti = NumCon(); + if (n > std::numeric_limits::max()) + return absl::InvalidArgumentError("arguments lb and ub too large"); + + MSK_appendcons(task.get(), (int)n); + std::vector bk(n); + for (ssize_t i = 0; i < n; ++i) { + bk[i] = (lb[i] > ub[i] + ? MSK_BK_RA + : (std::isfinite(lb[i]) + ? (std::isfinite(ub[i]) + ? (lb[i] < ub[i] ? MSK_BK_RA : MSK_BK_FX) + : MSK_BK_LO) + : (std::isfinite(ub[i]) ? MSK_BK_UP : MSK_BK_FR))); + } + + MSK_putconboundslice(task.get(), firsti, firsti + (int)n, bk.data(), + lb.data(), ub.data()); + return firsti; +} +absl::Status Mosek::PutVarType(VariableIndex j, bool is_integer) { + if (MSK_RES_OK != + MSK_putvartype(task.get(), j, + is_integer ? MSK_VAR_TYPE_INT : MSK_VAR_TYPE_CONT)) + return absl::InvalidArgumentError("Arguments j is invalid"); + return absl::OkStatus(); +} + +absl::Status Mosek::PutC(const std::vector& c) { + auto n = c.size(); + if (n > NumVar()) + return absl::InvalidArgumentError("Argument c is too large"); + for (int j = 0; j < n; ++j) MSK_putcj(task.get(), j, c[j]); + return absl::OkStatus(); +} + +void Mosek::PutCFix(double cfix) { MSK_putcfix(task.get(), cfix); } + +absl::Status Mosek::PutAIJList(const std::vector& subi, + const std::vector& subj, + const std::vector& valij) { + if (subi.size() != subj.size() || subi.size() != valij.size()) + return absl::InvalidArgumentError( + "Mismatching arguments subi, subj, valij"); + size_t n = subi.size(); + if (n > std::numeric_limits::max()) + return absl::InvalidArgumentError( + "Arguments subi, subj, valij are too long"); + + if (MSK_RES_OK != MSK_putaijlist(task.get(), (int)n, subi.data(), subj.data(), + valij.data())) + return absl::InvalidArgumentError("Invalid index argument subi or subj"); + return absl::OkStatus(); +} + +// Note: We implement indicator constraints as a disjunctive constraint of the +// form: [ indvar = (negate ? 1.0 : 0.0) ] +// OR +// [ indvar = (negate ? 0.0 : 1.0) +// lb <= Ax <= ub ] +// +absl::StatusOr +Mosek::AppendIndicatorConstraint(bool negate, VariableIndex indvar, + const std::vector& subj, + const std::vector& cof, double lb, + double ub) { + if (subj.size() != cof.size()) + return absl::InvalidArgumentError("Mismatching arguments subj, cof"); + + size_t n = subj.size(); + if (n > std::numeric_limits::max()) + return absl::InvalidArgumentError("Arguments subj or cof is too long"); + + int64_t ndjc, nafe; + MSK_getnumdjc(task.get(), &ndjc); + MSK_getnumafe(task.get(), &nafe); + + MSK_appendafes(task.get(), 2); + MSK_appenddjcs(task.get(), 1); + + MSK_putafefentry(task.get(), nafe, indvar, 1.0); + MSK_putafefrow(task.get(), nafe + 1, (int)n, subj.data(), cof.data()); + int64_t dom_eq, dom_lb, dom_ub; + MSK_appendrzerodomain(task.get(), 1, &dom_eq); + if (std::isfinite(lb)) { + MSK_appendrplusdomain(task.get(), 1, &dom_lb); + } else { + MSK_appendrdomain(task.get(), 1, &dom_lb); + } + if (std::isfinite(ub)) { + MSK_appendrminusdomain(task.get(), 1, &dom_ub); + } else { + MSK_appendrdomain(task.get(), 1, &dom_ub); + } + + int64_t afeidx[4] = {nafe, nafe, nafe + 1, nafe + 1}; + double b[4] = {negate ? 1.0 : 0.0, negate ? 0.0 : 1.0, lb, ub}; + int64_t domidxs[4] = {dom_eq, dom_eq, dom_lb, dom_ub}; + int64_t termsizes[2] = {1, 3}; + MSK_putdjc(task.get(), ndjc, 4, domidxs, 4, afeidx, b, 2, termsizes); + + return ndjc; +} +absl::Status Mosek::PutDJCName(DisjunctiveConstraintIndex djci, + const std::string& name) { + if (MSK_RES_OK != MSK_putdjcname(task.get(), djci, name.c_str())) + return absl::InvalidArgumentError("Invalid argument djci"); + return absl::OkStatus(); +} +absl::Status Mosek::PutACCName(ConeConstraintIndex acci, + const std::string& name) { + if (MSK_RES_OK != MSK_putaccname(task.get(), acci, name.c_str())) + return absl::InvalidArgumentError("Invalid argument acci"); + return absl::OkStatus(); +} + +absl::StatusOr Mosek::AppendConeConstraint( + ConeType ct, const std::vector& sizes, + const std::vector& subj, const std::vector& cof, + const std::vector& b) { + size_t n = sizes.size(); + size_t nnz = 0; + for (auto& nz : sizes) nnz += nz; + int64_t domidx; + + if (nnz != cof.size() || nnz != subj.size()) + return absl::InvalidArgumentError( + "Mismatching argument lengths of subj and cof"); + if (n != b.size()) + return absl::InvalidArgumentError("Mismatching argument lengths b and ptr"); + + switch (ct) { + case ConeType::SecondOrderCone: + MSK_appendquadraticconedomain(task.get(), n, &domidx); + break; + case ConeType::RotatedSecondOrderCone: + MSK_appendrquadraticconedomain(task.get(), n, &domidx); + break; + default: + return absl::InvalidArgumentError("Cone type not supported"); + } + + int64_t afei; + MSK_getnumafe(task.get(), &afei); + MSK_appendafes(task.get(), n); + + std::vector afeidxs(n); + for (ssize_t i = 0; i < n; ++i) afeidxs[i] = afei + i; + std::vector ptr(n); + ptr[0] = 0; + for (ssize_t i = 0; i < n; ++i) ptr[i + 1] = ptr[i] + sizes[i]; + + std::vector accb(n); + + int64_t acci; + MSK_getnumacc(task.get(), &acci); + MSK_appendaccseq(task.get(), domidx, n, afei, accb.data()); + MSK_putafefrowlist(task.get(), n, afeidxs.data(), sizes.data(), ptr.data(), + nnz, subj.data(), cof.data()); + for (ssize_t i = 0; i < n; ++i) MSK_putafeg(task.get(), afei + i, b[i]); + return acci; +} + +// Delete-ish +absl::Status Mosek::ClearVariable(VariableIndex j) { + if (MSK_RES_OK != MSK_putvarbound(task.get(), j, MSK_BK_FR, 0.0, 0.0)) + return absl::InvalidArgumentError("Invalid variable index j"); + return absl::OkStatus(); +} +absl::Status Mosek::ClearConstraint(ConstraintIndex i) { + if (MSK_RES_OK != MSK_putconbound(task.get(), i, MSK_BK_FR, 0.0, 0.0)) + return absl::InvalidArgumentError("Invalid constraint index i"); + int subj; + double cof; + MSK_putarow(task.get(), i, 0, &subj, &cof); + return absl::OkStatus(); +} +absl::Status Mosek::ClearConeConstraint(ConeConstraintIndex i) { + int64_t afeidxs; + double b; + if (MSK_RES_OK != MSK_putacc(task.get(), i, 0, 0, &afeidxs, &b)) + return absl::InvalidArgumentError("Invalid constraint index i"); + return absl::OkStatus(); +} +absl::Status Mosek::ClearDisjunctiveConstraint(DisjunctiveConstraintIndex i) { + int64_t i64s; + double f64s; + if (MSK_RES_OK != + MSK_putdjc(task.get(), i, 0, &i64s, 0, &i64s, &f64s, 0, &i64s)) + return absl::InvalidArgumentError("Invalid constraint index i"); + return absl::OkStatus(); +} + +// Update + +static MSKboundkeye merge_lower_bound(MSKboundkeye bk, double bl, double bu, + double b) { + switch (bk) { + case MSK_BK_FR: + return std::isfinite(b) ? MSK_BK_LO : MSK_BK_FR; + case MSK_BK_LO: + return std::isfinite(b) ? MSK_BK_LO : MSK_BK_FR; + case MSK_BK_UP: + return std::isfinite(b) ? (b < bu || b > bu ? MSK_BK_RA : MSK_BK_FX) + : MSK_BK_UP; + case MSK_BK_FX: + return std::isfinite(b) ? (b < bu || b > bu ? MSK_BK_RA : MSK_BK_FX) + : MSK_BK_UP; + case MSK_BK_RA: + return std::isfinite(b) ? (b < bu || b > bu ? MSK_BK_RA : MSK_BK_FX) + : MSK_BK_UP; + } + return MSK_BK_FX; +} + +static MSKboundkeye merge_upper_bound(MSKboundkeye bk, double bl, double bu, + double b) { + switch (bk) { + case MSK_BK_FR: + return std::isfinite(b) ? MSK_BK_UP : MSK_BK_FR; + case MSK_BK_UP: + return std::isfinite(b) ? MSK_BK_UP : MSK_BK_FR; + case MSK_BK_LO: + return std::isfinite(b) ? (b < bl || b > bl ? MSK_BK_RA : MSK_BK_FX) + : MSK_BK_UP; + case MSK_BK_FX: + return std::isfinite(b) ? (b < bl || b > bl ? MSK_BK_RA : MSK_BK_FX) + : MSK_BK_UP; + case MSK_BK_RA: + return std::isfinite(b) ? (b < bl || b > bl ? MSK_BK_RA : MSK_BK_FX) + : MSK_BK_UP; + } + return MSK_BK_FX; +} + +absl::Status Mosek::UpdateVariableLowerBound(VariableIndex j, double b) { + MSKboundkeye bk; + double bl, bu; + if (MSK_RES_OK != MSK_getvarbound(task.get(), j, &bk, &bl, &bu)) + return absl::InvalidArgumentError("Invalid variable index j"); + MSK_putvarbound(task.get(), j, merge_lower_bound(bk, bl, bu, b), b, bu); + return absl::OkStatus(); +} +absl::Status Mosek::UpdateVariableUpperBound(VariableIndex j, double b) { + MSKboundkeye bk; + double bl, bu; + if (MSK_RES_OK != MSK_getvarbound(task.get(), j, &bk, &bl, &bu)) + return absl::InvalidArgumentError("Invalid variable index j"); + MSK_putvarbound(task.get(), j, merge_upper_bound(bk, bl, bu, b), b, bu); + return absl::OkStatus(); +} +absl::Status Mosek::UpdateVariableType(VariableIndex j, bool is_integer) { + if (MSK_RES_OK != + MSK_putvartype(task.get(), j, + is_integer ? MSK_VAR_TYPE_INT : MSK_VAR_TYPE_CONT)) + return absl::InvalidArgumentError("Invalid variable index j"); + return absl::OkStatus(); +} +absl::Status Mosek::UpdateConstraintLowerBound(ConstraintIndex i, double b) { + MSKboundkeye bk; + double bl, bu; + if (MSK_RES_OK != MSK_getconbound(task.get(), i, &bk, &bl, &bu)) + return absl::InvalidArgumentError("Invalid constraint index i"); + MSK_putconbound(task.get(), i, merge_lower_bound(bk, bl, bu, b), b, bu); + return absl::OkStatus(); +} +absl::Status Mosek::UpdateConstraintUpperBound(ConstraintIndex i, double b) { + MSKboundkeye bk; + double bl, bu; + if (MSK_RES_OK != MSK_getconbound(task.get(), i, &bk, &bl, &bu)) + return absl::InvalidArgumentError("Invalid constraint index i"); + MSK_putconbound(task.get(), i, merge_upper_bound(bk, bl, bu, b), b, bu); + return absl::OkStatus(); +} +absl::Status Mosek::UpdateObjectiveSense(bool maximize) { + MSK_putobjsense(task.get(), maximize ? MSK_OBJECTIVE_SENSE_MAXIMIZE + : MSK_OBJECTIVE_SENSE_MINIMIZE); + return absl::OkStatus(); +} +absl::Status Mosek::UpdateObjective(double fixterm, + const std::vector& subj, + const std::vector& cof) { + if (subj.size() != cof.size()) + return absl::InvalidArgumentError( + "Mismatching argument lengths of subj and cof"); + if (subj.size() > std::numeric_limits::max()) + return absl::InvalidArgumentError("Argument subj and cof are too long"); + if (MSK_RES_OK != + MSK_putclist(task.get(), (int)subj.size(), subj.data(), cof.data())) + return absl::InvalidArgumentError("Invalid variable index in subj"); + MSK_putcfix(task.get(), fixterm); + return absl::OkStatus(); +} +absl::Status Mosek::UpdateA(const std::vector& subi, + const std::vector& subj, + const std::vector& cof) { + if (subi.size() != cof.size() || subj.size() != cof.size()) + return absl::InvalidArgumentError( + "Mismatching lengths of arguments subi, subj, and cof"); + if (MSK_RES_OK != MSK_putaijlist(task.get(), subi.size(), subi.data(), + subj.data(), cof.data())) + return absl::InvalidArgumentError( + "Invalid variable or constraint index in subi or subj"); + return absl::OkStatus(); +} + +void Mosek::WriteData(const char* filename) const { + MSK_writedata(task.get(), filename); +} + +absl::StatusOr Mosek::Optimize() { + MSKrescodee trm; + MSKrescodee r = MSK_optimizetrm(task.get(), &trm); + if (MSK_RES_OK != r) return absl::InternalError("Optimization failed"); + + return trm; +} + +void Mosek::message_callback(MSKuserhandle_t handle, const char* msg) { + MessageCallback& msg_cb = *reinterpret_cast(handle); + if (msg_cb) (msg_cb)(msg); +} +int Mosek::info_callback(MSKtask_t task, MSKuserhandle_t h, + MSKcallbackcodee code, const double* dinf, + const int* iinf, const int64_t* liinf) { + auto& info_cb = *reinterpret_cast(h); + + bool interrupt = false; + if (info_cb) + interrupt = info_cb(code, absl::Span(dinf, MSK_DINF_END), + absl::Span(iinf, MSK_IINF_END), + absl::Span(liinf, MSK_LIINF_END)); + return interrupt ? 1 : 0; +} + +absl::StatusOr Mosek::Optimize(MessageCallback msg_cb, + InfoCallback info_cb) { + MSKrescodee trm; + + auto _cleanup = absl::MakeCleanup([&]() { + MSK_linkfunctotaskstream(task.get(), MSK_STREAM_LOG, nullptr, nullptr); + MSK_putcallbackfunc(task.get(), nullptr, nullptr); + }); + + bool interrupt = false; + + if (info_cb) MSK_putcallbackfunc(task.get(), info_callback, &info_cb); + if (msg_cb) + MSK_linkfunctotaskstream(task.get(), MSK_STREAM_LOG, &msg_cb, + message_callback); + + MSKrescodee r = MSK_optimizetrm(task.get(), &trm); + if (MSK_RES_OK != r) return absl::InternalError("Optimization failed"); + + return trm; +} + +std::tuple Mosek::LastError() const { + int64_t msglen; + MSKrescodee r; + if (MSK_RES_OK != MSK_getlasterror64(task.get(), &r, 0, &msglen, nullptr)) { + return {"", "", MSK_RES_ERR_UNKNOWN}; + } else { + std::vector msg(msglen + 1); + char buf[MSK_MAX_STR_LEN]; + MSK_getlasterror64(task.get(), &r, msglen, &msglen, msg.data()); + msg[msglen] = 0; + + MSK_rescodetostr(r, buf); + return {msg.data(), buf, r}; + } +} + +double Mosek::GetPrimalObj(MSKsoltypee whichsol) const { + if (!SolutionDef(whichsol)) return 0.0; + double val; + MSK_getprimalobj(task.get(), whichsol, &val); + return val; +} +double Mosek::GetDualObj(MSKsoltypee whichsol) const { + if (!SolutionDef(whichsol)) return 0.0; + double val; + MSK_getdualobj(task.get(), whichsol, &val); + return val; +} +void Mosek::GetXX(MSKsoltypee whichsol, std::vector& xx) const { + xx.clear(); + if (!SolutionDef(whichsol)) return; + xx.resize(NumVar()); + MSK_getxx(task.get(), whichsol, xx.data()); +} +void Mosek::GetSLX(MSKsoltypee whichsol, std::vector& slx) const { + slx.clear(); + if (!SolutionDef(whichsol)) return; + slx.resize(NumVar()); + MSK_getslx(task.get(), whichsol, slx.data()); +} +void Mosek::GetSUX(MSKsoltypee whichsol, std::vector& sux) const { + sux.clear(); + if (!SolutionDef(whichsol)) return; + sux.resize(NumVar()); + MSK_getsux(task.get(), whichsol, sux.data()); +} +void Mosek::GetSLC(MSKsoltypee whichsol, std::vector& slc) const { + slc.clear(); + if (!SolutionDef(whichsol)) return; + slc.resize(NumVar()); + MSK_getslc(task.get(), whichsol, slc.data()); +} +void Mosek::GetSUC(MSKsoltypee whichsol, std::vector& suc) const { + suc.clear(); + if (!SolutionDef(whichsol)) return; + suc.resize(NumVar()); + MSK_getsuc(task.get(), whichsol, suc.data()); +} +void Mosek::GetY(MSKsoltypee whichsol, std::vector& y) const { + y.clear(); + if (!SolutionDef(whichsol)) return; + y.resize(NumVar()); + MSK_gety(task.get(), whichsol, y.data()); +} +void Mosek::GetSKX(MSKsoltypee whichsol, std::vector& skx) const { + skx.clear(); + if (!SolutionDef(whichsol)) return; + skx.resize(NumVar()); + MSK_getskx(task.get(), whichsol, skx.data()); +} +void Mosek::GetSKC(MSKsoltypee whichsol, std::vector& skc) const { + skc.clear(); + if (!SolutionDef(whichsol)) return; + skc.resize(NumVar()); + MSK_getskc(task.get(), whichsol, skc.data()); +} + +void Mosek::PutParam(MSKiparame par, int value) { + MSK_putintparam(task.get(), par, value); +} +void Mosek::PutParam(MSKdparame par, double value) { + MSK_putdouparam(task.get(), par, value); +} +// Query +int Mosek::NumVar() const { + int n; + MSK_getnumvar(task.get(), &n); + return n; +} +int Mosek::NumCon() const { + int n; + MSK_getnumcon(task.get(), &n); + return n; +} +bool Mosek::IsMaximize() const { + MSKobjsensee sense; + MSK_getobjsense(task.get(), &sense); + return sense == MSK_OBJECTIVE_SENSE_MAXIMIZE; +} +double Mosek::GetParam(MSKdparame dpar) const { + double parval = 0.0; + MSK_getdouparam(task.get(), dpar, &parval); + return parval; +} +int Mosek::GetParam(MSKiparame ipar) const { + int parval; + MSK_getintparam(task.get(), ipar, &parval); + return parval; +} + +void Mosek::delete_msk_task_func(MSKtask_t t) { MSK_deletetask(&t); } +} // namespace operations_research::math_opt diff --git a/ortools/math_opt/solvers/mosek/mosekwrp.h b/ortools/math_opt/solvers/mosek/mosekwrp.h new file mode 100644 index 00000000000..cae3f5041e5 --- /dev/null +++ b/ortools/math_opt/solvers/mosek/mosekwrp.h @@ -0,0 +1,180 @@ +#ifndef OR_TOOLS_MATH_OPT_SOLVERS_MOSEK_MOSEKWRP_H_ +#define OR_TOOLS_MATH_OPT_SOLVERS_MOSEK_MOSEKWRP_H_ + +#include +#include +#include + +#include "absl/container/flat_hash_map.h" +#include "absl/status/statusor.h" +#include "absl/types/span.h" +#include "mosek.h" + +namespace operations_research::math_opt { + +namespace mosek { +enum class SolSta : int { + UNKNOWN = MSK_SOL_STA_UNKNOWN, + OPTIMAL = MSK_SOL_STA_OPTIMAL, + PRIM_FEAS = MSK_SOL_STA_PRIM_FEAS, + DUAL_FEAS = MSK_SOL_STA_DUAL_FEAS, + PRIM_AND_DUAL_FEAS = MSK_SOL_STA_PRIM_AND_DUAL_FEAS, + PRIM_INFEAS_CER = MSK_SOL_STA_PRIM_INFEAS_CER, + DUAL_INFEAS_CER = MSK_SOL_STA_DUAL_INFEAS_CER, + PRIM_ILLPOSED_CER = MSK_SOL_STA_PRIM_ILLPOSED_CER, + DUAL_ILLPOSED_CER = MSK_SOL_STA_DUAL_ILLPOSED_CER, + INTEGER_OPTIMAL = MSK_SOL_STA_INTEGER_OPTIMAL +}; +enum class ProSta : int { + UNKNOWN = MSK_PRO_STA_UNKNOWN, + PRIM_AND_DUAL_FEAS = MSK_PRO_STA_PRIM_AND_DUAL_FEAS, + PRIM_FEAS = MSK_PRO_STA_PRIM_FEAS, + DUAL_FEAS = MSK_PRO_STA_DUAL_FEAS, + PRIM_INFEAS = MSK_PRO_STA_PRIM_INFEAS, + DUAL_INFEAS = MSK_PRO_STA_DUAL_INFEAS, + PRIM_AND_DUAL_INFEAS = MSK_PRO_STA_PRIM_AND_DUAL_INFEAS, + ILL_POSED = MSK_PRO_STA_ILL_POSED, + PRIM_INFEAS_OR_UNBOUNDED = MSK_PRO_STA_PRIM_INFEAS_OR_UNBOUNDED +}; + +std::ostream& operator<<(std::ostream& s, SolSta solsta); +std::ostream& operator<<(std::ostream& s, ProSta prosta); +} // namespace mosek +class Mosek { + public: + enum class ConeType { SecondOrderCone, RotatedSecondOrderCone }; + + using MessageCallback = std::function; + using InfoCallback = + std::function, + absl::Span, absl::Span)>; + + typedef int32_t VariableIndex; + typedef int32_t ConstraintIndex; + typedef int64_t DisjunctiveConstraintIndex; + typedef int64_t ConeConstraintIndex; + + Mosek(Mosek&& m); + Mosek(Mosek& m) = delete; + + static Mosek* Create(); + + void PutName(const std::string& name); + void PutObjName(const std::string& name); + void PutVarName(VariableIndex j, const std::string& name); + void PutConName(ConstraintIndex j, const std::string& name); + void PutObjectiveSense(bool maximize); + + absl::StatusOr AppendVars(const std::vector& lb, + const std::vector& ub); + absl::StatusOr AppendCons(const std::vector& lb, + const std::vector& ub); + absl::Status PutVarType(VariableIndex j, bool is_integer); + + absl::Status PutC(const std::vector& c); + void PutCFix(double cfix); + + absl::Status PutAIJList(const std::vector& subi, + const std::vector& subj, + const std::vector& valij); + + absl::StatusOr AppendIndicatorConstraint( + bool negate, VariableIndex indvar, const std::vector& subj, + const std::vector& cof, double lb, double ub); + absl::Status PutDJCName(DisjunctiveConstraintIndex djci, + const std::string& name); + absl::Status PutACCName(ConeConstraintIndex acci, const std::string& name); + + absl::StatusOr AppendConeConstraint( + ConeType ct, const std::vector& sizes, + const std::vector& subj, const std::vector& cof, + const std::vector& b); + + // Delete-ish + absl::Status ClearVariable(VariableIndex j); + absl::Status ClearConstraint(ConstraintIndex i); + absl::Status ClearConeConstraint(ConeConstraintIndex i); + absl::Status ClearDisjunctiveConstraint(DisjunctiveConstraintIndex i); + + // Update + + absl::Status UpdateVariableLowerBound(VariableIndex j, double b); + absl::Status UpdateVariableUpperBound(VariableIndex j, double b); + absl::Status UpdateVariableType(VariableIndex j, bool is_integer); + absl::Status UpdateConstraintLowerBound(ConstraintIndex i, double b); + absl::Status UpdateConstraintUpperBound(ConstraintIndex i, double b); + absl::Status UpdateObjectiveSense(bool maximize); + absl::Status UpdateObjective(double fixterm, + const std::vector& subj, + const std::vector& cof); + absl::Status UpdateA(const std::vector& subi, + const std::vector& subj, + const std::vector& cof); + + // Solve + absl::StatusOr Optimize(); + + absl::StatusOr Optimize(MessageCallback msg_cb, + InfoCallback info_cb); + + // Parameters + + void PutParam(MSKiparame par, int value); + void PutParam(MSKdparame par, double value); + + // Query + int NumVar() const; + int NumCon() const; + bool IsMaximize() const; + double GetParam(MSKdparame dpar) const; + int GetParam(MSKiparame ipar) const; + + bool SolutionDef(MSKsoltypee which) const { + MSKbooleant soldef; + MSK_solutiondef(task.get(), which, &soldef); + return soldef != 0; + } + mosek::ProSta GetProSta(MSKsoltypee which) const { + MSKprostae prosta; + MSK_getprosta(task.get(), which, &prosta); + return (mosek::ProSta)prosta; + } + mosek::SolSta GetSolSta(MSKsoltypee which) const { + MSKsolstae solsta; + MSK_getsolsta(task.get(), which, &solsta); + return (mosek::SolSta)solsta; + } + + std::tuple LastError() const; + + double GetPrimalObj(MSKsoltypee whichsol) const; + double GetDualObj(MSKsoltypee whichsol) const; + + void GetXX(MSKsoltypee whichsol, std::vector& xx) const; + void GetSLX(MSKsoltypee whichsol, std::vector& slx) const; + void GetSUX(MSKsoltypee whichsol, std::vector& sux) const; + void GetSLC(MSKsoltypee whichsol, std::vector& slc) const; + void GetSUC(MSKsoltypee whichsol, std::vector& suc) const; + void GetY(MSKsoltypee whichsol, std::vector& y) const; + void GetSKX(MSKsoltypee whichsol, std::vector& skx) const; + void GetSKC(MSKsoltypee whichsol, std::vector& skc) const; + + // Other + void WriteData(const char* filename) const; + + private: + static void delete_msk_task_func(MSKtask_t); + + Mosek(MSKtask_t task); + + std::unique_ptr task; + + static void message_callback(MSKuserhandle_t handle, const char* msg); + static int info_callback(MSKtask_t task, MSKuserhandle_t h, + MSKcallbackcodee code, const double* dinf, + const int* iinf, const int64_t* liinf); +}; + +} // namespace operations_research::math_opt + +#endif diff --git a/ortools/math_opt/solvers/mosek_solver.cc b/ortools/math_opt/solvers/mosek_solver.cc new file mode 100644 index 00000000000..60c29f8012d --- /dev/null +++ b/ortools/math_opt/solvers/mosek_solver.cc @@ -0,0 +1,1181 @@ +// Copyright 2010-2024 Google LLskip_xx_zeros C +// 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 "mosek_solver.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// #include "absl/algorithm/container.h" +// #include "absl/cleanup/cleanup.h" +// #include "absl/container/flat_hash_map.h" +#include "absl/log/check.h" +// #include "absl/memory/memory.h" +#include "absl/cleanup/cleanup.h" +#include "absl/status/status.h" +#include "absl/status/statusor.h" +#include "absl/strings/string_view.h" +// #include "absl/time/clock.h" +// #include "absl/time/time.h" +#include "absl/types/span.h" +#include "ortools/base/protoutil.h" +#include "ortools/base/status_builder.h" +#include "ortools/base/status_macros.h" +#include "ortools/math_opt/callback.pb.h" +#include "ortools/math_opt/core/empty_bounds.h" +#include "ortools/math_opt/core/inverted_bounds.h" +#include "ortools/math_opt/core/math_opt_proto_utils.h" +#include "ortools/math_opt/core/solver_interface.h" +#include "ortools/math_opt/core/sorted.h" +#include "ortools/math_opt/core/sparse_vector_view.h" +#include "ortools/math_opt/infeasible_subsystem.pb.h" +#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/message_callback_data.h" +#include "ortools/math_opt/solvers/mosek.pb.h" +#include "ortools/util/solve_interrupter.h" +#include "ortools/util/status_macros.h" + +namespace operations_research::math_opt { +namespace { + +constexpr SupportedProblemStructures kMosekSupportedStructures = { + .integer_variables = SupportType::kSupported, + .second_order_cone_constraints = SupportType::kSupported, + .indicator_constraints = SupportType::kSupported, +}; + +} // namespace + +std::ostream& mosek::operator<<(std::ostream& s, SolSta solsta) { + switch (solsta) { + case SolSta::UNKNOWN: + s << "UNKNOWN"; + break; + case SolSta::OPTIMAL: + s << "OPTIMAL"; + break; + case SolSta::PRIM_FEAS: + s << "PRIM_FEAS"; + break; + case SolSta::DUAL_FEAS: + s << "DUAL_FEAS"; + break; + case SolSta::PRIM_AND_DUAL_FEAS: + s << "PRIM_AND_DUAL_FEAS"; + break; + case SolSta::PRIM_INFEAS_CER: + s << "PRIM_INFEAS_CER"; + break; + case SolSta::DUAL_INFEAS_CER: + s << "DUAL_INFEAS_CER"; + break; + case SolSta::PRIM_ILLPOSED_CER: + s << "PRIM_ILLPOSED_CER"; + break; + case SolSta::DUAL_ILLPOSED_CER: + s << "DUAL_ILLPOSED_CER"; + break; + case SolSta::INTEGER_OPTIMAL: + s << "INTEGER_OPTIMAL"; + break; + } + return s; +} + +std::ostream& mosek::operator<<(std::ostream& s, ProSta prosta) { + switch (prosta) { + case ProSta::UNKNOWN: + s << "UNKNOWN"; + break; + case ProSta::PRIM_AND_DUAL_FEAS: + s << "PRIM_AND_DUAL_FEAS"; + break; + case ProSta::PRIM_FEAS: + s << "PRIM_FEAS"; + break; + case ProSta::DUAL_FEAS: + s << "DUAL_FEAS"; + break; + case ProSta::PRIM_INFEAS: + s << "PRIM_INFEAS"; + break; + case ProSta::DUAL_INFEAS: + s << "DUAL_INFEAS"; + break; + case ProSta::PRIM_AND_DUAL_INFEAS: + s << "PRIM_AND_DUAL_INFEAS"; + break; + case ProSta::ILL_POSED: + s << "ILL_POSED"; + break; + case ProSta::PRIM_INFEAS_OR_UNBOUNDED: + s << "PRIM_INFEAS_OR_UNBOUNDED"; + break; + } + return s; +} + +absl::Status MosekSolver::AddVariables(const VariablesProto& vars) { + int num_vars = vars.ids_size(); + int firstvar = msk.NumVar(); + + { + int i = 0; + for (const auto& v : vars.ids()) { + variable_map[v] = firstvar + i; + ++i; + } + } + + std::vector lbx(num_vars); + std::vector ubx(num_vars); + { + int i = 0; + for (const auto& v : vars.lower_bounds()) lbx[i++] = v; + } + { + int i = 0; + for (const auto& v : vars.upper_bounds()) ubx[i++] = v; + } + + { + auto r = msk.AppendVars(lbx, ubx); + if (!r.ok()) return r.status(); + } + + { + int i = 0; + for (const bool is_integer : vars.integers()) { + if (is_integer) { + RETURN_IF_ERROR(msk.PutVarType(variable_map[i], true)); + } + ++i; + } + } + { + int i = 0; + for (const auto& name : vars.names()) { + msk.PutVarName(firstvar + i, name.c_str()); + ++i; + } + } + return absl::OkStatus(); +} // MosekSolver::AddVariables + +absl::Status MosekSolver::ReplaceObjective(const ObjectiveProto& obj) { + msk.PutObjName(obj.name()); + RETURN_IF_ERROR(msk.UpdateObjectiveSense(obj.maximize())); + auto objcof = obj.linear_coefficients(); + msk.PutCFix(obj.offset()); + auto num_vars = msk.NumVar(); + std::vector c(num_vars); + auto n = objcof.ids_size(); + for (int64_t i = 0; i < n; ++i) { + c[objcof.ids(i)] = objcof.values(i); + } + RETURN_IF_ERROR(msk.PutC(c)); + return absl::OkStatus(); +} // MosekSolver::ReplaceObjective + +absl::Status MosekSolver::AddConstraints(const LinearConstraintsProto& cons, + const SparseDoubleMatrixProto& adata) { + int firstcon = msk.NumCon(); + auto numcon = cons.ids_size(); + { + int i = 0; + for (const auto& id : cons.ids()) { + linconstr_map[id] = i; + ++i; + } + } + std::vector clb(numcon); + std::vector cub(numcon); + { + int i = 0; + for (const auto& b : cons.lower_bounds()) clb[i++] = b; + } + { + int i = 0; + for (const auto& b : cons.upper_bounds()) cub[i++] = b; + } + { + auto r = msk.AppendCons(clb, cub); + if (!r.ok()) return r.status(); + } + + { + int i = 0; + for (const auto& name : cons.names()) { + msk.PutConName(firstcon + i, name.c_str()); + ++i; + } + } + + size_t nnz = adata.row_ids_size(); + std::vector subj; + subj.reserve(nnz); + std::vector subi; + subi.reserve(nnz); + std::vector valij; + valij.reserve(nnz); + + for (const auto id : adata.row_ids()) subi.push_back(linconstr_map[id]); + for (const auto id : adata.column_ids()) subj.push_back(variable_map[id]); + for (const auto c : adata.coefficients()) valij.push_back(c); + RETURN_IF_ERROR(msk.PutAIJList(subi, subj, valij)); + + return absl::OkStatus(); +} +absl::Status MosekSolver::AddConstraints(const LinearConstraintsProto& cons) { + int firstcon = msk.NumCon(); + auto numcon = cons.ids_size(); + { + int i = 0; + for (const auto& id : cons.ids()) { + linconstr_map[id] = i; + ++i; + } + } + std::vector clb(numcon); + std::vector cub(numcon); + { + int i = 0; + for (const auto& b : cons.lower_bounds()) clb[i++] = b; + } + { + int i = 0; + for (const auto& b : cons.upper_bounds()) cub[i++] = b; + } + { + auto r = msk.AppendCons(clb, cub); + if (!r.ok()) return r.status(); + } + + { + int i = 0; + for (const auto& name : cons.names()) { + msk.PutConName(firstcon + i, name.c_str()); + ++i; + } + } + + return absl::OkStatus(); +} // MosekSolver::AddConstraints + +absl::Status MosekSolver::AddIndicatorConstraints( + const ::google::protobuf::Map& cons) { + int i = 0; + std::vector subj; + std::vector cof; + for (const auto& [id, con] : cons) { + indconstr_map[id] = i++; + Mosek::VariableIndex indvar = indconstr_map[con.indicator_id()]; + + subj.clear(); + subj.reserve(con.expression().ids_size()); + cof.clear(); + cof.reserve(con.expression().ids_size()); + + for (auto id : con.expression().ids()) { + subj.push_back(variable_map[id]); + } + for (auto c : con.expression().values()) { + cof.push_back(c); + } + + auto djci = + msk.AppendIndicatorConstraint(con.activate_on_zero(), indvar, subj, cof, + con.lower_bound(), con.upper_bound()); + if (!djci.ok()) { + return djci.status(); + } + + RETURN_IF_ERROR(msk.PutDJCName(*djci, con.name())); + } + return absl::OkStatus(); +} // MosekSolver::AddIndicatorConstraints + +absl::Status MosekSolver::AddConicConstraints( + const ::google::protobuf::Map& + cons) { + std::vector subj; + std::vector cof; + std::vector sizes; + std::vector b; + + sizes.reserve(cons.size()); + for (const auto& [idx, con] : cons) { + auto& expr0 = con.upper_bound(); + int64_t totalnnz = expr0.ids_size(); + for (const auto& lexp : con.arguments_to_norm()) { + totalnnz += lexp.ids_size(); + } + + subj.reserve(totalnnz); + cof.reserve(totalnnz); + b.push_back(expr0.offset()); + + for (const auto& id : expr0.ids()) { + subj.push_back(variable_map[id]); + } + for (auto c : expr0.coefficients()) { + cof.push_back(c); + } + + for (const auto& expri : con.arguments_to_norm()) { + sizes.push_back(expri.ids_size()); + for (const auto& id : expri.ids()) { + subj.push_back(variable_map[id]); + } + for (auto c : expri.coefficients()) { + cof.push_back(c); + } + b.push_back(expri.offset()); + } + + auto acci = msk.AppendConeConstraint(Mosek::ConeType::SecondOrderCone, + sizes, subj, cof, b); + if (!acci.ok()) { + return acci.status(); + } + + RETURN_IF_ERROR(msk.PutACCName(*acci, con.name())); + } + return absl::OkStatus(); +} + +absl::StatusOr MosekSolver::Update(const ModelUpdateProto& model_update) { + for (auto id : model_update.deleted_variable_ids()) { + variable_map.erase(id); + RETURN_IF_ERROR(msk.ClearVariable(variable_map[id])); + } + for (auto id : model_update.deleted_linear_constraint_ids()) { + linconstr_map.erase(id); + RETURN_IF_ERROR(msk.ClearConstraint(linconstr_map[id])); + } + for (auto id : model_update.second_order_cone_constraint_updates() + .deleted_constraint_ids()) { + coneconstr_map.erase(id); + RETURN_IF_ERROR(msk.ClearConeConstraint(coneconstr_map[id])); + } + for (auto id : + model_update.indicator_constraint_updates().deleted_constraint_ids()) { + indconstr_map.erase(id); + RETURN_IF_ERROR(msk.ClearDisjunctiveConstraint(indconstr_map[id])); + } + + RETURN_IF_ERROR(AddVariables(model_update.new_variables())); + RETURN_IF_ERROR(UpdateVariables(model_update.variable_updates())); + RETURN_IF_ERROR(AddConstraints(model_update.new_linear_constraints())); + RETURN_IF_ERROR( + UpdateConstraints(model_update.linear_constraint_updates(), + model_update.linear_constraint_matrix_updates())); + + RETURN_IF_ERROR(UpdateObjective(model_update.objective_updates())); + RETURN_IF_ERROR(AddConicConstraints( + model_update.second_order_cone_constraint_updates().new_constraints())); + RETURN_IF_ERROR(AddIndicatorConstraints( + model_update.indicator_constraint_updates().new_constraints())); + // RETURN_IF_ERROR(UpdateIndicatorConstraint(conupd)); + return true; +} + +absl::Status MosekSolver::UpdateVariables(const VariableUpdatesProto& varupds) { + // std::cout << "MosekSolver::UpdateVariables()" << std::endl; + for (int64_t i = 0, n = varupds.lower_bounds().ids_size(); i < n; ++i) { + RETURN_IF_ERROR(msk.UpdateVariableLowerBound( + variable_map[varupds.lower_bounds().ids(i)], + varupds.lower_bounds().values(i))); + } + for (int64_t i = 0, n = varupds.upper_bounds().ids_size(); i < n; ++i) { + RETURN_IF_ERROR(msk.UpdateVariableUpperBound( + variable_map[varupds.upper_bounds().ids(i)], + varupds.upper_bounds().values(i))); + } + for (int64_t i = 0, n = varupds.integers().ids_size(); i < n; ++i) { + RETURN_IF_ERROR( + msk.UpdateVariableType(variable_map[varupds.upper_bounds().ids(i)], + varupds.integers().values(i))); + } + return absl::OkStatus(); +} +absl::Status MosekSolver::UpdateConstraints( + const LinearConstraintUpdatesProto& conupds, + const SparseDoubleMatrixProto& lincofupds) { + for (int64_t i = 0, n = conupds.lower_bounds().ids_size(); i < n; ++i) { + RETURN_IF_ERROR(msk.UpdateConstraintLowerBound( + linconstr_map[conupds.lower_bounds().ids(i)], + conupds.lower_bounds().values(i))); + } + for (int64_t i = 0, n = conupds.upper_bounds().ids_size(); i < n; ++i) { + RETURN_IF_ERROR(msk.UpdateConstraintUpperBound( + linconstr_map[conupds.upper_bounds().ids(i)], + conupds.upper_bounds().values(i))); + } + + size_t n = lincofupds.row_ids_size(); + std::vector subi(n); + std::vector subj(n); + std::vector valij(lincofupds.coefficients().begin(), + lincofupds.coefficients().end()); + { + int i = 0; + for (auto id : lincofupds.row_ids()) { + subi[i] = linconstr_map[id]; + ++i; + } + } + { + int i = 0; + for (auto id : lincofupds.column_ids()) { + subj[i] = variable_map[id]; + ++i; + } + } + + RETURN_IF_ERROR(msk.UpdateA(subi, subj, valij)); + return absl::OkStatus(); +} +absl::Status MosekSolver::UpdateObjective( + const ObjectiveUpdatesProto& objupds) { + const auto& vals = objupds.linear_coefficients(); + std::vector cof(vals.values().begin(), vals.values().end()); + std::vector subj; + subj.reserve(cof.size()); + for (auto id : objupds.linear_coefficients().ids()) + subj.push_back(variable_map[id]); + + RETURN_IF_ERROR(msk.UpdateObjectiveSense(objupds.direction_update())); + RETURN_IF_ERROR(msk.UpdateObjective(objupds.offset_update(), subj, cof)); + + return absl::OkStatus(); +} +absl::Status MosekSolver::UpdateConstraint( + const SecondOrderConeConstraintUpdatesProto& conupds) { + for (auto id : conupds.deleted_constraint_ids()) { + RETURN_IF_ERROR(msk.ClearConeConstraint(coneconstr_map[id])); + } + + RETURN_IF_ERROR(AddConicConstraints(conupds.new_constraints())); + + return absl::OkStatus(); +} + +absl::Status MosekSolver::UpdateConstraint( + const IndicatorConstraintUpdatesProto& conupds) { + for (auto id : conupds.deleted_constraint_ids()) { + RETURN_IF_ERROR(msk.ClearDisjunctiveConstraint(indconstr_map[id])); + } + + RETURN_IF_ERROR(AddIndicatorConstraints(conupds.new_constraints())); + + return absl::OkStatus(); +} + +absl::StatusOr> MosekSolver::New( + const ModelProto& model, const InitArgs&) { + RETURN_IF_ERROR(ModelIsSupported(model, kMosekSupportedStructures, "Mosek")); + + if (!model.auxiliary_objectives().empty()) + return util::InvalidArgumentErrorBuilder() + << "Mosek does not support multi-objective models"; + if (!model.objective().quadratic_coefficients().row_ids().empty()) { + return util::InvalidArgumentErrorBuilder() + << "Mosek does not support models with quadratic objectives"; + } + if (!model.quadratic_constraints().empty()) { + return util::InvalidArgumentErrorBuilder() + << "Mosek does not support models with quadratic constraints"; + } + if (!model.sos1_constraints().empty() || !model.sos2_constraints().empty()) { + return util::InvalidArgumentErrorBuilder() + << "Mosek does not support models with SOS constraints"; + } + + std::unique_ptr msk(Mosek::Create()); + std::unique_ptr mskslv(new MosekSolver(std::move(*msk))); + mskslv->msk.PutName(model.name()); + + RETURN_IF_ERROR(mskslv->AddVariables(model.variables())); + RETURN_IF_ERROR(mskslv->ReplaceObjective(model.objective())); + RETURN_IF_ERROR(mskslv->AddConstraints(model.linear_constraints(), + model.linear_constraint_matrix())); + RETURN_IF_ERROR( + mskslv->AddIndicatorConstraints(model.indicator_constraints())); + + std::unique_ptr res(std::move(mskslv)); + + return res; +} + +MosekSolver::MosekSolver(Mosek&& msk) : msk(std::move(msk)) {} + +absl::StatusOr MosekSolver::PrimalSolution( + MSKsoltypee whichsol, const std::vector& ordered_var_ids, + bool skip_zero_values) { + auto solsta = msk.GetSolSta(whichsol); + PrimalSolutionProto sol; + switch (solsta) { + case mosek::SolSta::OPTIMAL: + case mosek::SolSta::INTEGER_OPTIMAL: + case mosek::SolSta::PRIM_AND_DUAL_FEAS: + case mosek::SolSta::PRIM_FEAS: + sol.set_feasibility_status(SolutionStatusProto::SOLUTION_STATUS_FEASIBLE); + { + sol.set_objective_value(msk.GetPrimalObj(whichsol)); + std::vector xx; + msk.GetXX(whichsol, xx); + SparseDoubleVectorProto vals; + + for (auto k : ordered_var_ids) { + auto v = xx[variable_map[k]]; + if (!skip_zero_values || v < 0 || v > 0) { + vals.add_ids(k); + vals.add_values(v); + } + } + + *sol.mutable_variable_values() = std::move(vals); + } + break; + default: + return absl::NotFoundError("Primal solution not available"); + } + return std::move(sol); +} +absl::StatusOr MosekSolver::DualSolution( + MSKsoltypee whichsol, const std::vector& ordered_y_ids, + bool skip_y_zeros, const std::vector& ordered_yx_ids, + bool skip_yx_zeros) { + auto solsta = msk.GetSolSta(whichsol); + DualSolutionProto sol; + switch (solsta) { + case mosek::SolSta::OPTIMAL: + case mosek::SolSta::PRIM_AND_DUAL_FEAS: + case mosek::SolSta::DUAL_FEAS: { + sol.set_objective_value(msk.GetPrimalObj(whichsol)); + sol.set_feasibility_status(SolutionStatusProto::SOLUTION_STATUS_FEASIBLE); + std::vector keys; + keys.reserve(std::max(variable_map.size(), linconstr_map.size())); + { + std::vector slx; + msk.GetSLX(whichsol, slx); + std::vector sux; + msk.GetSUX(whichsol, sux); + SparseDoubleVectorProto vals; + + for (auto k : ordered_yx_ids) { + auto j = variable_map[k]; + auto v = slx[j] - sux[j]; + if (!skip_yx_zeros || v < 0.0 || v > 0.0) { + vals.add_ids(k); + vals.add_values(v); + } + } + *sol.mutable_reduced_costs() = std::move(vals); + } + { + std::vector y; + msk.GetY(whichsol, y); + SparseDoubleVectorProto vals; + for (auto k : ordered_y_ids) { + auto v = y[linconstr_map[k]]; + if (!skip_y_zeros || v < 0 || v > 0) { + vals.add_ids(k); + vals.add_values(v); + } + } + + *sol.mutable_dual_values() = std::move(vals); + } + } break; + default: + return absl::NotFoundError("Primal solution not available"); + } + return std::move(sol); +} +absl::StatusOr MosekSolver::Solution( + MSKsoltypee whichsol, const std::vector& ordered_xc_ids, + const std::vector& ordered_xx_ids, bool skip_xx_zeros, + const std::vector& ordered_y_ids, bool skip_y_zeros, + const std::vector& ordered_yx_ids, bool skip_yx_zeros) { + // std::cout << "MosekSolver::Solution()" << std::endl; + SolutionProto sol; + { + auto r = PrimalSolution(whichsol, ordered_xx_ids, skip_xx_zeros); + if (r.ok()) *sol.mutable_primal_solution() = std::move(*r); + } + { + auto r = DualSolution(whichsol, ordered_y_ids, skip_y_zeros, ordered_yx_ids, + skip_yx_zeros); + if (r.ok()) *sol.mutable_dual_solution() = std::move(*r); + } + + if (whichsol == MSK_SOL_BAS) { + // std::cout << "MosekSolver::Solution(): Basis!" << std::endl; + BasisProto bas; + SparseBasisStatusVector csta; + SparseBasisStatusVector xsta; + std::vector sk; + msk.GetSKX(whichsol, sk); + + for (auto k : ordered_xx_ids) { + auto v = variable_map[k]; + xsta.add_ids(k); + switch (sk[v]) { + case MSK_SK_LOW: + xsta.add_values(BasisStatusProto::BASIS_STATUS_AT_LOWER_BOUND); + break; + case MSK_SK_UPR: + xsta.add_values(BasisStatusProto::BASIS_STATUS_AT_UPPER_BOUND); + break; + case MSK_SK_FIX: + xsta.add_values(BasisStatusProto::BASIS_STATUS_FIXED_VALUE); + break; + case MSK_SK_BAS: + xsta.add_values(BasisStatusProto::BASIS_STATUS_BASIC); + break; + case MSK_SK_INF: + case MSK_SK_SUPBAS: + case MSK_SK_UNK: + xsta.add_values(BasisStatusProto::BASIS_STATUS_UNSPECIFIED); + break; + } + } + sk.clear(); + msk.GetSKC(whichsol, sk); + for (auto k : ordered_xc_ids) { + auto v = linconstr_map[k]; + csta.add_ids(k); + switch (sk[v]) { + case MSK_SK_LOW: + csta.add_values(BasisStatusProto::BASIS_STATUS_AT_LOWER_BOUND); + break; + case MSK_SK_UPR: + csta.add_values(BasisStatusProto::BASIS_STATUS_AT_UPPER_BOUND); + break; + case MSK_SK_FIX: + csta.add_values(BasisStatusProto::BASIS_STATUS_FIXED_VALUE); + break; + case MSK_SK_BAS: + csta.add_values(BasisStatusProto::BASIS_STATUS_BASIC); + break; + case MSK_SK_INF: + case MSK_SK_SUPBAS: + case MSK_SK_UNK: + csta.add_values(BasisStatusProto::BASIS_STATUS_UNSPECIFIED); + break; + } + } + *bas.mutable_variable_status() = std::move(xsta); + *bas.mutable_constraint_status() = std::move(csta); + + auto solsta = msk.GetSolSta(whichsol); + switch (solsta) { + case mosek::SolSta::OPTIMAL: + case mosek::SolSta::INTEGER_OPTIMAL: + case mosek::SolSta::PRIM_AND_DUAL_FEAS: + case mosek::SolSta::PRIM_FEAS: + bas.set_basic_dual_feasibility( + SolutionStatusProto::SOLUTION_STATUS_FEASIBLE); + break; + default: + bas.set_basic_dual_feasibility( + SolutionStatusProto::SOLUTION_STATUS_UNSPECIFIED); + break; + } + + *sol.mutable_basis() = std::move(bas); + } + return std::move(sol); +} + +absl::StatusOr MosekSolver::PrimalRay( + MSKsoltypee whichsol, const std::vector& ordered_xx_ids, + bool skip_xx_zeros) { + auto solsta = msk.GetSolSta(whichsol); + if (solsta == mosek::SolSta::DUAL_INFEAS_CER) + return absl::NotFoundError("Certificate not available"); + + std::vector xx; + msk.GetXX(whichsol, xx); + PrimalRayProto ray; + SparseDoubleVectorProto data; + for (auto k : ordered_xx_ids) { + auto v = xx[variable_map[k]]; + if (!skip_xx_zeros || v < 0 || v > 0) { + data.add_ids(k); + data.add_values(v); + } + } + *ray.mutable_variable_values() = data; + return ray; +} + +absl::StatusOr MosekSolver::DualRay( + MSKsoltypee whichsol, const std::vector& ordered_y_ids, + bool skip_y_zeros, const std::vector& ordered_yx_ids, + bool skip_yx_zeros) { + auto solsta = msk.GetSolSta(whichsol); + + if (solsta == mosek::SolSta::PRIM_INFEAS_CER) + return absl::NotFoundError("Certificate not available"); + + std::vector slx; + msk.GetSLX(whichsol, slx); + std::vector sux; + msk.GetSUX(whichsol, slx); + std::vector y; + msk.GetY(whichsol, y); + DualRayProto ray; + SparseDoubleVectorProto xdata; + SparseDoubleVectorProto cdata; + for (auto k : ordered_yx_ids) { + auto j = variable_map[k]; + auto v = slx[j] - sux[j]; + if (!skip_yx_zeros || v < 0 || v > 0) { + xdata.add_ids(k); + xdata.add_values(v); + } + } + for (auto& k : ordered_y_ids) { + auto v = y[linconstr_map[k]]; + if (!skip_y_zeros || v < 0 || v > 0) { + cdata.add_ids(k); + cdata.add_values(v); + } + } + *ray.mutable_dual_values() = xdata; + *ray.mutable_reduced_costs() = cdata; + return ray; +} + +absl::StatusOr MosekSolver::Solve( + const SolveParametersProto& parameters, // solver settings + const ModelSolveParametersProto& model_parameters, + MessageCallback message_cb, + const CallbackRegistrationProto& callback_registration, + Callback cb, // using Callback = std::function, from base_solver.h + const SolveInterrupter* const solve_interrupter) { + // Solve parameters that we support: + // - google.protobuf.Duration time_limit + // - optional int64 iteration_limit + // - optional int64 node_limit + // - optional double cutoff_limit + // - bool enable_output + // - optional int32 threads + // - optional double absolute_gap_tolerance + // - optional double relative_gap_tolerance + // - LPAlgorithmProto lp_algorithm + // Solve parameters that we may support: + // - optional double best_bound_limit + // - optional double objective_limit + // Solve parameters that we do not support: + // - optional int32 solution_pool_size + // - optional int32 solution_limit + // - optional int32 random_seed + // - EmphasisProto presolve + // - EmphasisProto cuts + // - EmphasisProto heuristics + // - EmphasisProto scaling + + double dpar_optimizer_max_time = msk.GetParam(MSK_DPAR_OPTIMIZER_MAX_TIME); + int ipar_intpnt_max_iterations = msk.GetParam(MSK_IPAR_INTPNT_MAX_ITERATIONS); + int ipar_sim_max_iterations = msk.GetParam(MSK_IPAR_SIM_MAX_ITERATIONS); + double dpar_upper_obj_cut = msk.GetParam(MSK_DPAR_UPPER_OBJ_CUT); + double dpar_lower_obj_cut = msk.GetParam(MSK_DPAR_LOWER_OBJ_CUT); + int ipar_num_threads = msk.GetParam(MSK_IPAR_NUM_THREADS); + double dpar_mio_tol_abs_gap = msk.GetParam(MSK_DPAR_MIO_TOL_ABS_GAP); + double dpar_mio_tol_rel_gap = msk.GetParam(MSK_DPAR_MIO_TOL_REL_GAP); + double dpar_intpnt_tol_rel_gap = msk.GetParam(MSK_DPAR_INTPNT_TOL_REL_GAP); + double dpar_intpnt_co_tol_rel_gap = + msk.GetParam(MSK_DPAR_INTPNT_CO_TOL_REL_GAP); + int ipar_optimizer = msk.GetParam(MSK_IPAR_OPTIMIZER); + + auto _guard_reset_params = absl::MakeCleanup([&]() { + msk.PutParam(MSK_DPAR_OPTIMIZER_MAX_TIME, dpar_optimizer_max_time); + msk.PutParam(MSK_IPAR_INTPNT_MAX_ITERATIONS, ipar_intpnt_max_iterations); + msk.PutParam(MSK_IPAR_SIM_MAX_ITERATIONS, ipar_sim_max_iterations); + msk.PutParam(MSK_DPAR_UPPER_OBJ_CUT, dpar_upper_obj_cut); + msk.PutParam(MSK_DPAR_LOWER_OBJ_CUT, dpar_lower_obj_cut); + msk.PutParam(MSK_IPAR_NUM_THREADS, ipar_num_threads); + msk.PutParam(MSK_DPAR_MIO_TOL_ABS_GAP, dpar_mio_tol_abs_gap); + msk.PutParam(MSK_DPAR_MIO_TOL_REL_GAP, dpar_mio_tol_rel_gap); + msk.PutParam(MSK_DPAR_INTPNT_TOL_REL_GAP, dpar_intpnt_tol_rel_gap); + msk.PutParam(MSK_DPAR_INTPNT_CO_TOL_REL_GAP, dpar_intpnt_co_tol_rel_gap); + }); + + if (parameters.has_time_limit()) { + OR_ASSIGN_OR_RETURN3( + const absl::Duration time_limit, + util_time::DecodeGoogleApiProto(parameters.time_limit()), + _ << "invalid time_limit value for HiGHS."); + msk.PutParam(MSK_DPAR_OPTIMIZER_MAX_TIME, + absl::ToDoubleSeconds(time_limit)); + } + + if (parameters.has_iteration_limit()) { + const int iter_limit = parameters.iteration_limit(); + + msk.PutParam(MSK_IPAR_INTPNT_MAX_ITERATIONS, iter_limit); + msk.PutParam(MSK_IPAR_SIM_MAX_ITERATIONS, iter_limit); + } + + // Not supported in MOSEK 10.2 + // int ipar_mio_ + // if (parameters.has_node_limit()) { + // ASSIGN_OR_RETURN( + // const int node_limit, + // SafeIntCast(parameters.node_limit(), "node_limit")); + // msk.PutIntParam(MSK_IPAR_MIO__MAX_NODES, node_limit); + //} + + // Not supported by MOSEK? + // if (parameters.has_cutoff_limit()) { + //} + if (parameters.has_objective_limit()) { + if (msk.IsMaximize()) + msk.PutParam(MSK_DPAR_UPPER_OBJ_CUT, parameters.cutoff_limit()); + else + msk.PutParam(MSK_DPAR_LOWER_OBJ_CUT, parameters.cutoff_limit()); + } + + if (parameters.has_threads()) { + msk.PutParam(MSK_IPAR_NUM_THREADS, parameters.threads()); + } + + if (parameters.has_absolute_gap_tolerance()) { + msk.PutParam(MSK_DPAR_MIO_TOL_ABS_GAP, parameters.absolute_gap_tolerance()); + } + + if (parameters.has_relative_gap_tolerance()) { + msk.PutParam(MSK_DPAR_INTPNT_TOL_REL_GAP, + parameters.absolute_gap_tolerance()); + msk.PutParam(MSK_DPAR_INTPNT_CO_TOL_REL_GAP, + parameters.absolute_gap_tolerance()); + msk.PutParam(MSK_DPAR_MIO_TOL_REL_GAP, parameters.absolute_gap_tolerance()); + } + + switch (parameters.lp_algorithm()) { + case LP_ALGORITHM_BARRIER: + msk.PutParam(MSK_IPAR_OPTIMIZER, MSK_OPTIMIZER_INTPNT); + break; + case LP_ALGORITHM_DUAL_SIMPLEX: + msk.PutParam(MSK_IPAR_OPTIMIZER, MSK_OPTIMIZER_DUAL_SIMPLEX); + break; + case LP_ALGORITHM_PRIMAL_SIMPLEX: + msk.PutParam(MSK_IPAR_OPTIMIZER, MSK_OPTIMIZER_PRIMAL_SIMPLEX); + break; + default: + // use default auto select, usually intpnt + msk.PutParam(MSK_IPAR_OPTIMIZER, MSK_OPTIMIZER_FREE); + break; + } + + // TODO: parameter enable_output + + bool skip_xx_zeros = + model_parameters.variable_values_filter().skip_zero_values(); + bool skip_y_zeros = model_parameters.dual_values_filter().skip_zero_values(); + bool skip_yx_zeros = + model_parameters.reduced_costs_filter().skip_zero_values(); + bool filter_ids = model_parameters.variable_values_filter().filter_by_ids(); + + std::vector ordered_xc_ids; + std::vector ordered_xx_ids; + std::vector ordered_y_ids; + std::vector ordered_yx_ids; + + ordered_xc_ids.reserve(linconstr_map.size()); + for (auto [id, idx] : linconstr_map) ordered_xc_ids.push_back(id); + std::sort(ordered_xc_ids.begin(), ordered_xc_ids.end()); + + if (!skip_xx_zeros) { + ordered_xx_ids.reserve(variable_map.size()); + for (auto [id, idx] : variable_map) { + ordered_xx_ids.push_back(id); + } + std::sort(ordered_xx_ids.begin(), ordered_xx_ids.end()); + } else { + ordered_xx_ids.reserve( + model_parameters.variable_values_filter().filtered_ids().size()); + for (auto id : model_parameters.variable_values_filter().filtered_ids()) { + if (variable_map.contains(id)) ordered_xx_ids.push_back(id); + } + } + + if (!model_parameters.dual_values_filter().filter_by_ids()) { + ordered_y_ids.reserve(linconstr_map.size()); + for (auto [id, idx] : linconstr_map) ordered_y_ids.push_back(id); + std::sort(ordered_y_ids.begin(), ordered_y_ids.end()); + } else { + ordered_y_ids.reserve( + model_parameters.dual_values_filter().filtered_ids().size()); + for (auto id : model_parameters.dual_values_filter().filtered_ids()) + ordered_y_ids.push_back(id); + } + + if (!model_parameters.reduced_costs_filter().filter_by_ids()) { + ordered_yx_ids.reserve(linconstr_map.size()); + for (auto [id, idx] : variable_map) ordered_yx_ids.push_back(id); + std::sort(ordered_yx_ids.begin(), ordered_yx_ids.end()); + } else { + ordered_yx_ids.reserve( + model_parameters.reduced_costs_filter().filtered_ids().size()); + for (auto id : model_parameters.reduced_costs_filter().filtered_ids()) + ordered_yx_ids.push_back(id); + } + + MSKrescodee trm; + { + BufferedMessageCallback bmsg_cb(message_cb); + // TODO: Use model_parameters + auto r = msk.Optimize( + [&](const std::string& msg) { bmsg_cb.OnMessage(msg); }, + [&](MSKcallbackcodee code, absl::Span dinf, + absl::Span iinf, absl::Span liinf) { + if (cb) { + CallbackDataProto cbdata; + switch (code) { + case MSK_CALLBACK_IM_SIMPLEX: + cbdata.mutable_simplex_stats()->set_iteration_count( + liinf[MSK_LIINF_SIMPLEX_ITER]); + cbdata.mutable_simplex_stats()->set_objective_value( + dinf[MSK_DINF_SIM_OBJ]); + // cbdata.mutable_simplex_stats(/->set_primal_infeasibility(...); + // cbdata.mutable_simplex_stats()->set_dual_infeasibility(...); + // cbdata.mutable_simplex_stats()->is_perturbed(...); + cbdata.set_event(CALLBACK_EVENT_SIMPLEX); + break; + case MSK_CALLBACK_IM_MIO: + cbdata.mutable_mip_stats()->set_primal_bound( + dinf[MSK_DINF_MIO_OBJ_BOUND]); + // cbdata.mutable_mip_stats()->set_dual_bound(...); + cbdata.mutable_mip_stats()->set_explored_nodes( + iinf[MSK_IINF_MIO_NUM_SOLVED_NODES]); + // cbdata.mutable_mip_stats()->set_open_nodes(...); + cbdata.mutable_mip_stats()->set_simplex_iterations( + liinf[MSK_LIINF_MIO_SIMPLEX_ITER]); + // cbdata.mutable_mip_stats()->set_number_of_solutions_found(...); + // cbdata.mutable_mip_stats()->set_cutting_planes_in_lp(...); + cbdata.set_event(CALLBACK_EVENT_MIP); + break; + case MSK_CALLBACK_NEW_INT_MIO: + cbdata.set_event(CALLBACK_EVENT_MIP_SOLUTION); + { + std::vector xx; + msk.GetXX(MSK_SOL_ITG, xx); + + SparseDoubleVectorProto primal; + + for (auto id : ordered_xx_ids) { + auto v = xx[variable_map[id]]; + if (!skip_xx_zeros || v > 0.0 || v < 0.0) { + primal.add_ids(id); + primal.add_values(v); + } + } + *cbdata.mutable_primal_solution_vector() = primal; + } + break; + case MSK_CALLBACK_IM_PRESOLVE: + cbdata.set_event(CALLBACK_EVENT_PRESOLVE); + break; + case MSK_CALLBACK_IM_CONIC: + case MSK_CALLBACK_IM_INTPNT: + cbdata.mutable_barrier_stats()->set_iteration_count( + liinf[MSK_IINF_INTPNT_ITER]); + cbdata.mutable_barrier_stats()->set_primal_objective( + dinf[MSK_DINF_INTPNT_PRIMAL_OBJ]); + cbdata.mutable_barrier_stats()->set_dual_objective( + dinf[MSK_DINF_INTPNT_DUAL_OBJ]); + // cbdata.mutable_barrier_stats()->set_complementarity(...); + // cbdata.mutable_barrier_stats()->set_primal_infeasibility(...); + // cbdata.mutable_barrier_stats()->set_dual_infeasibility(...); + + cbdata.set_event(CALLBACK_EVENT_BARRIER); + break; + default: + cbdata.set_event(CALLBACK_EVENT_UNSPECIFIED); + break; + } + + auto r = cb(cbdata); + if (r.ok()) { + return r->terminate(); + } + } + return false; + }); + // msk.WriteData("__test.ptf"); + // std::cout << "MosekSolver::Solve() optimize -> " << r << std::endl; + if (!r.ok()) return r.status(); + trm = *r; + } + + MSKsoltypee whichsol{}; + bool soldef = true; + if (msk.SolutionDef(MSK_SOL_ITG)) { + whichsol = MSK_SOL_ITG; + } else if (msk.SolutionDef(MSK_SOL_BAS)) { + whichsol = MSK_SOL_BAS; + } else if (msk.SolutionDef(MSK_SOL_ITR)) { + whichsol = MSK_SOL_ITR; + } else { + soldef = false; + } + + TerminationProto trmp; + mosek::ProSta prosta{}; + mosek::SolSta solsta{}; + if (!soldef) { + auto [msg, name, code] = msk.LastError(); + trmp = TerminateForReason( + msk.IsMaximize(), + TerminationReasonProto::TERMINATION_REASON_NO_SOLUTION_FOUND, msg); + } else { + // std::cout << "MosekSolver::Solve() solution is defined " << whichsol << + // std::endl; + prosta = msk.GetProSta(whichsol); + solsta = msk.GetSolSta(whichsol); + + // Attempt to determine TerminationProto from Mosek Termination code, + // problem status and solution status. + + if (solsta == mosek::SolSta::OPTIMAL || + solsta == mosek::SolSta::INTEGER_OPTIMAL) { + // std::cout << "MosekSolver::Solve() trmp = Optimal! " << std::endl; + trmp = OptimalTerminationProto(msk.GetPrimalObj(whichsol), + msk.GetDualObj(whichsol), ""); + } else if (solsta == mosek::SolSta::PRIM_INFEAS_CER) + trmp = InfeasibleTerminationProto( + msk.IsMaximize(), + FeasibilityStatusProto::FEASIBILITY_STATUS_FEASIBLE); + else if (prosta == mosek::ProSta::PRIM_INFEAS_OR_UNBOUNDED) + trmp = InfeasibleOrUnboundedTerminationProto(msk.IsMaximize()); + else if (solsta == mosek::SolSta::DUAL_INFEAS_CER) + trmp = UnboundedTerminationProto(msk.IsMaximize()); + else if (solsta == mosek::SolSta::PRIM_AND_DUAL_FEAS || + solsta == mosek::SolSta::PRIM_FEAS) { + LimitProto lim = LimitProto::LIMIT_UNSPECIFIED; + switch (trm) { + case MSK_RES_TRM_MAX_ITERATIONS: + lim = LimitProto::LIMIT_ITERATION; + break; + case MSK_RES_TRM_MAX_TIME: + lim = LimitProto::LIMIT_TIME; + break; + case MSK_RES_TRM_NUM_MAX_NUM_INT_SOLUTIONS: + lim = LimitProto::LIMIT_SOLUTION; + break; +#if MSK_VERSION_MAJOR >= 11 + case MSK_RES_TRM_SERVER_MAX_MEMORY: + lim = LimitProto::LIMIT_MEMORY; + break; +#endif + // LIMIT_CUTOFF + case MSK_RES_TRM_OBJECTIVE_RANGE: + lim = LimitProto::LIMIT_OBJECTIVE; + break; + case MSK_RES_TRM_NUMERICAL_PROBLEM: + lim = LimitProto::LIMIT_NORM; + break; + case MSK_RES_TRM_USER_CALLBACK: + lim = LimitProto::LIMIT_INTERRUPTED; + break; + case MSK_RES_TRM_STALL: + lim = LimitProto::LIMIT_SLOW_PROGRESS; + break; + default: + lim = LimitProto::LIMIT_OTHER; + break; + } + if (solsta == mosek::SolSta::PRIM_AND_DUAL_FEAS) + trmp = FeasibleTerminationProto(msk.IsMaximize(), lim, + msk.GetPrimalObj(whichsol), + msk.GetDualObj(whichsol)); + else + trmp = FeasibleTerminationProto( + msk.IsMaximize(), lim, msk.GetPrimalObj(whichsol), std::nullopt); + } else { + trmp = NoSolutionFoundTerminationProto(msk.IsMaximize(), + LimitProto::LIMIT_UNSPECIFIED); + } + } + + SolveResultProto result; + *result.mutable_termination() = trmp; + + if (soldef) { + // TODO: Use model_parameters + switch (solsta) { + case mosek::SolSta::OPTIMAL: + case mosek::SolSta::INTEGER_OPTIMAL: + case mosek::SolSta::PRIM_FEAS: + case mosek::SolSta::DUAL_FEAS: + case mosek::SolSta::PRIM_AND_DUAL_FEAS: { + auto r = Solution(whichsol, ordered_xc_ids, ordered_xx_ids, + skip_xx_zeros, ordered_y_ids, skip_y_zeros, + ordered_yx_ids, skip_yx_zeros); + if (r.ok()) { + *result.add_solutions() = std::move(*r); + } + } break; + case mosek::SolSta::DUAL_INFEAS_CER: { + auto r = PrimalRay(whichsol, ordered_xx_ids, skip_xx_zeros); + if (r.ok()) { + *result.add_primal_rays() = std::move(*r); + } + } break; + case mosek::SolSta::PRIM_INFEAS_CER: { + auto r = DualRay(whichsol, ordered_y_ids, skip_y_zeros, ordered_yx_ids, + skip_yx_zeros); + if (r.ok()) { + *result.add_dual_rays() = std::move(*r); + } + } break; + case mosek::SolSta::PRIM_ILLPOSED_CER: + case mosek::SolSta::DUAL_ILLPOSED_CER: + case mosek::SolSta::UNKNOWN: + break; + } + } + return result; +} + +absl::StatusOr +MosekSolver::ComputeInfeasibleSubsystem(const SolveParametersProto&, + MessageCallback, + const SolveInterrupter*) { + return absl::UnimplementedError( + "MOSEK does not yet support computing an infeasible subsystem"); +} + +MATH_OPT_REGISTER_SOLVER(SOLVER_TYPE_MOSEK, MosekSolver::New); + +} // namespace operations_research::math_opt diff --git a/ortools/math_opt/solvers/mosek_solver.h b/ortools/math_opt/solvers/mosek_solver.h new file mode 100644 index 00000000000..12aeb992169 --- /dev/null +++ b/ortools/math_opt/solvers/mosek_solver.h @@ -0,0 +1,188 @@ +// Copyright 2010-2024 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_MOSEK_SOLVER_H_ +#define OR_TOOLS_MATH_OPT_SOLVERS_MOSEK_SOLVER_H_ + +#include +#include +#include +#include +#include + +#include "absl/container/flat_hash_map.h" +#include "absl/status/status.h" +#include "absl/status/statusor.h" +#include "mosek/mosekwrp.h" +#include "ortools/base/status_builder.h" +#include "ortools/math_opt/callback.pb.h" +#include "ortools/math_opt/core/inverted_bounds.h" +#include "ortools/math_opt/core/solver_interface.h" +#include "ortools/math_opt/model.pb.h" +#include "ortools/math_opt/model_parameters.pb.h" +#include "ortools/math_opt/model_update.pb.h" +#include "ortools/math_opt/parameters.pb.h" +#include "ortools/math_opt/result.pb.h" +#include "ortools/math_opt/solution.pb.h" +#include "ortools/util/solve_interrupter.h" + +namespace operations_research::math_opt { + +class MosekSolver : public SolverInterface { + public: + static absl::StatusOr> New( + const ModelProto& model, const InitArgs& init_args); + + absl::StatusOr Solve( + const SolveParametersProto& parameters, + const ModelSolveParametersProto& model_parameters, + MessageCallback message_cb, + const CallbackRegistrationProto& callback_registration, Callback cb, + const SolveInterrupter* interrupter) override; + absl::StatusOr Update(const ModelUpdateProto& model_update) override; + absl::StatusOr + ComputeInfeasibleSubsystem(const SolveParametersProto& parameters, + MessageCallback message_cb, + const SolveInterrupter* interrupter) override; + + private: + Mosek msk; + absl::flat_hash_map variable_map; + absl::flat_hash_map linconstr_map; + absl::flat_hash_map coneconstr_map; + absl::flat_hash_map indconstr_map; + + absl::Status ReplaceObjective(const ObjectiveProto& obj); + absl::Status AddVariables(const VariablesProto& vars); + absl::Status AddConstraints(const LinearConstraintsProto& cons); + absl::Status AddConstraints(const LinearConstraintsProto& cons, + const SparseDoubleMatrixProto& lincofupds); + absl::Status AddIndicatorConstraints( + const ::google::protobuf::Map& cons); + absl::Status AddConicConstraints( + const ::google::protobuf::Map& + cons); + + absl::Status UpdateVariables(const VariableUpdatesProto& varupds); + absl::Status UpdateConstraints(const LinearConstraintUpdatesProto& conupds, + const SparseDoubleMatrixProto& lincofupds); + absl::Status UpdateObjective(const ObjectiveUpdatesProto& objupds); + absl::Status UpdateConstraint( + const SecondOrderConeConstraintUpdatesProto& conupds); + absl::Status UpdateConstraint(const IndicatorConstraintUpdatesProto& conupds); + + absl::StatusOr PrimalSolution( + MSKsoltypee whichsol, const std::vector& ordered_xx_ids, + bool skip_zeros); + absl::StatusOr DualSolution( + MSKsoltypee whichsol, const std::vector& ordered_y_ids, + bool skip_y_zeros, const std::vector& ordered_yx_ids, + bool skip_yx_zeros); + absl::StatusOr Solution( + MSKsoltypee whichsol, const std::vector& ordered_xc_ids, + const std::vector& ordered_xx_ids, bool skip_xx_zeros, + const std::vector& ordered_y_ids, bool skip_y_zeros, + const std::vector& ordered_yx_ids, bool skip_yx_zeros); + absl::StatusOr PrimalRay( + MSKsoltypee whichsol, const std::vector& ordered_xx_ids, + bool skip_zeros); + absl::StatusOr DualRay( + MSKsoltypee whichsol, const std::vector& ordered_y_ids, + bool skip_y_zeros, const std::vector& ordered_yx_ids, + bool skip_yx_zeros); + + MosekSolver(Mosek&& msk); + MosekSolver(MosekSolver&) = delete; + +#if 0 + struct SolutionClaims { + bool mosek_returned_primal_feasible_solution = false; + bool mosek_returned_dual_feasible_solution = false; + bool mosek_returned_primal_ray = false; + bool mosek_returned_dual_ray = false; + }; + struct SolutionsAndClaims { + std::vector solutions; + SolutionClaims solution_claims; + }; + + + absl::StatusOr PrimalRayReturned() const; + absl::StatusOr DualRayReturned() const; + + // Will append solutions and rays to result_out. No other modifications are + // made to result_out. Requires that mosek_->getInfo is validated. + absl::StatusOr ExtractSolutionAndRays( + const ModelSolveParametersProto& model_params); + + static absl::StatusOr PrimalFeasibilityStatus( + SolutionClaims solution_claims); + static absl::StatusOr DualFeasibilityStatus( + const MosekInfo& mosek_info, bool is_integer, + SolutionClaims solution_claims); + static absl::StatusOr MakeTermination( + MosekModelStatus mosek_model_status, const MosekInfo& mosek_info, + bool is_integer, bool had_node_limit, bool had_solution_limit, + bool is_maximize, SolutionClaims solution_claims); + + // Returns the current basis if it is available and MathOpt can represent it + // (all kNonBasic values can be made more precise, see b/272767311). + absl::StatusOr> ExtractBasis(); + + template + absl::Status EnsureOneEntryPerVariable(const std::vector& vec) { + if (vec.size() != variable_data_.size()) { + return util::InvalidArgumentErrorBuilder() + << "expected one entry per variable, but model had " + << variable_data_.size() << " variables and found " << vec.size() + << " elements"; + } + return absl::OkStatus(); + } + + template + absl::Status EnsureOneEntryPerLinearConstraint(const std::vector& vec) { + if (vec.size() != lin_con_data_.size()) { + return util::InvalidArgumentErrorBuilder() + << "expected one entry per linear constraint, but model had " + << lin_con_data_.size() << " linear constraints and found " + << vec.size() << " elements"; + } + return absl::OkStatus(); + } + + // Returns a SolveResult for when HiGHS returns the model status + // MosekModelStatus::kModelEmpty. This happens on models that have no + // variables, but may still have (potentially infeasible) linear constraints + // and an objective offset. + // + // Assumes that there are no inverted linear constraint bounds. + static SolveResultProto ResultForMosekModelStatusModelEmpty( + bool is_maximize, double objective_offset, + const absl::flat_hash_map& lin_con_data); + + InvertedBounds ListInvertedBounds(); + + std::unique_ptr mosek_; + + // Key is the mathopt id, value.index is the variable index in HiGHS. + absl::flat_hash_map variable_data_; + + // Key is the mathopt id, value.index is the linear constraint index in HiGHS. + absl::flat_hash_map lin_con_data_; +#endif +}; + +} // namespace operations_research::math_opt + +#endif // OR_TOOLS_MATH_OPT_SOLVERS_MOSEK_SOLVER_H_ diff --git a/ortools/service/v1/mathopt/parameters.proto b/ortools/service/v1/mathopt/parameters.proto index 5dc0e4e524c..4fb30b7fd8e 100644 --- a/ortools/service/v1/mathopt/parameters.proto +++ b/ortools/service/v1/mathopt/parameters.proto @@ -99,6 +99,11 @@ enum SolverTypeProto { // Slow/not recommended for production. Not an LP solver (no dual information // returned). SOLVER_TYPE_SANTORINI = 11; + + // Moske solver (third party). + // + // Supports LP, MIP, and conic quadratic. Requires a license. + SOLVER_TYPE_MOSEK = 12; } // Selects an algorithm for solving linear programs. From f29855e37f85fd1458e19807c5cb6a37dcaf9664 Mon Sep 17 00:00:00 2001 From: ulfw Date: Mon, 28 Oct 2024 12:16:14 +0100 Subject: [PATCH 02/10] fix: mosek link isuse --- cmake/check_deps.cmake | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/cmake/check_deps.cmake b/cmake/check_deps.cmake index a0b4ba90b86..43b3833d2e0 100644 --- a/cmake/check_deps.cmake +++ b/cmake/check_deps.cmake @@ -117,8 +117,11 @@ if(USE_CPLEX AND NOT TARGET CPLEX::CPLEX) endif() # Check optional Dependencies -if(USE_MOSEK AND NOT TARGET mosek::mosek) - message(FATAL_ERROR "Target mosek::mosek not available.") +if(USE_MOSEK) + if (NOT TARGET mosek::mosek) + message(FATAL_ERROR "Target mosek::mosek not available.") + endif() + set(MOSEK_DEPS mosek::mosek) endif() # CXX Test From 96b7cc37fdee0e77cec0cff911b8959f218428aa Mon Sep 17 00:00:00 2001 From: ulfw Date: Mon, 28 Oct 2024 15:15:12 +0100 Subject: [PATCH 03/10] add support for quadratic terms in MOSEK --- ortools/math_opt/solvers/mosek/mosekwrp.cc | 86 ++++++++- ortools/math_opt/solvers/mosek/mosekwrp.h | 15 ++ ortools/math_opt/solvers/mosek_solver.cc | 203 ++++++++++++++++++++- ortools/math_opt/solvers/mosek_solver.h | 85 +-------- 4 files changed, 297 insertions(+), 92 deletions(-) diff --git a/ortools/math_opt/solvers/mosek/mosekwrp.cc b/ortools/math_opt/solvers/mosek/mosekwrp.cc index f1f64420a06..ecb4a8f553d 100644 --- a/ortools/math_opt/solvers/mosek/mosekwrp.cc +++ b/ortools/math_opt/solvers/mosek/mosekwrp.cc @@ -70,16 +70,27 @@ absl::StatusOr Mosek::AppendVars( lb.data(), ub.data()); return firstj; } +absl::StatusOr Mosek::AppendCons( + double lb, double ub) { + return append_cons(1,&lb,&ub); +} + absl::StatusOr Mosek::AppendCons( const std::vector& lb, const std::vector& ub) { if (lb.size() != ub.size()) return absl::InvalidArgumentError("Mismatching lengths of lb and ub"); - size_t n = lb.size(); + int n = (int)lb.size(); + if (n != lb.size()) + return absl::InvalidArgumentError("Number of added constraints is too large"); + return append_cons(n,lb.data(),ub.data()); +} + +absl::StatusOr Mosek::append_cons( + int n, const double * lb, const double * ub) { int firsti = NumCon(); - if (n > std::numeric_limits::max()) - return absl::InvalidArgumentError("arguments lb and ub too large"); - MSK_appendcons(task.get(), (int)n); + if (MSK_RES_OK != MSK_appendcons(task.get(), n)) + return absl::InternalError("Maximum number of internal constraints exceeded"); std::vector bk(n); for (ssize_t i = 0; i < n; ++i) { bk[i] = (lb[i] > ub[i] @@ -91,10 +102,10 @@ absl::StatusOr Mosek::AppendCons( : (std::isfinite(ub[i]) ? MSK_BK_UP : MSK_BK_FR))); } - MSK_putconboundslice(task.get(), firsti, firsti + (int)n, bk.data(), - lb.data(), ub.data()); + MSK_putconboundslice(task.get(), firsti, firsti + (int)n, bk.data(), lb, ub); return firsti; } + absl::Status Mosek::PutVarType(VariableIndex j, bool is_integer) { if (MSK_RES_OK != MSK_putvartype(task.get(), j, @@ -103,6 +114,43 @@ absl::Status Mosek::PutVarType(VariableIndex j, bool is_integer) { return absl::OkStatus(); } +absl::Status Mosek::PutQObj(const std::vector & subk, + const std::vector & subl, + const std::vector & valkl) { + + if (subk.size() != subl.size() || + subk.size() != valkl.size()) + return absl::InvalidArgumentError("Mismatching argument lengths of subk, subl, valkl"); + + int qnnz = (int)subk.size(); + if (qnnz != subk.size()) + return absl::InvalidArgumentError("Number of quadratic objective nonzeros is too large"); + + if (MSK_RES_OK != + MSK_putqobj(task.get(),qnnz,subk.data(),subl.data(),valkl.data())) + return absl::InternalError("Error setting quadratic objective terms"); + return absl::OkStatus(); +} +absl::Status Mosek::UpdateQObjEntries(const std::vector & subk, + const std::vector & subl, + const std::vector & valkl) { + MSKrescodee r = MSK_RES_OK; + if (subk.size() != subl.size() || + subk.size() != valkl.size()) + return absl::InvalidArgumentError("Mismatching argument lengths"); + for (auto [kit,lit,cit] = std::make_tuple(subk.begin(), subl.begin(),valkl.begin()); + MSK_RES_OK == r && + kit != subk.end() && + lit != subl.end() && + cit != valkl.end(); + ++kit, ++lit, ++cit) + r = MSK_putqobjij(task.get(),*kit,*lit,*cit); + + if (r != MSK_RES_OK) + return absl::InternalError("Failed to input quadratic objective nonzeros"); + return absl::OkStatus(); +} + absl::Status Mosek::PutC(const std::vector& c) { auto n = c.size(); if (n > NumVar()) @@ -113,6 +161,31 @@ absl::Status Mosek::PutC(const std::vector& c) { void Mosek::PutCFix(double cfix) { MSK_putcfix(task.get(), cfix); } +absl::Status Mosek::PutQCon(int i, const std::vector& subk, + const std::vector& subl, + const std::vector& cof) { + if (subk.size() != subl.size() || + subk.size() != cof.size()) + return absl::InvalidArgumentError("Mismatching argument sizes"); + if (MSK_RES_OK != + MSK_putqconk(task.get(),i,subk.size(),subk.data(),subl.data(),cof.data())) + return absl::InternalError("Failed to set quadratic nonzeros in constraint"); + + return absl::OkStatus(); +} + +absl::Status Mosek::PutARow(int i, const std::vector& subj, + const std::vector& cof) { + if (subj.size() != cof.size()) + return absl::InvalidArgumentError("Mismatching argument lengths"); + if (subj.size() > std::numeric_limits::max()) + return absl::InvalidArgumentError("Number of nonzeros added too large"); + int nnz = (int)subj.size(); + + if (MSK_RES_OK != MSK_putarow(task.get(),i,nnz,subj.data(),cof.data())) + return absl::InternalError("Failed to add row nonzeros"); + return absl::OkStatus(); +} absl::Status Mosek::PutAIJList(const std::vector& subi, const std::vector& subj, const std::vector& valij) { @@ -250,6 +323,7 @@ absl::Status Mosek::ClearConstraint(ConstraintIndex i) { int subj; double cof; MSK_putarow(task.get(), i, 0, &subj, &cof); + MSK_putqconk(task.get(), i, 0, &subj,&subj,&cof); return absl::OkStatus(); } absl::Status Mosek::ClearConeConstraint(ConeConstraintIndex i) { diff --git a/ortools/math_opt/solvers/mosek/mosekwrp.h b/ortools/math_opt/solvers/mosek/mosekwrp.h index cae3f5041e5..dbde1eb4bb8 100644 --- a/ortools/math_opt/solvers/mosek/mosekwrp.h +++ b/ortools/math_opt/solvers/mosek/mosekwrp.h @@ -67,13 +67,25 @@ class Mosek { absl::StatusOr AppendVars(const std::vector& lb, const std::vector& ub); + absl::StatusOr AppendCons(double lb, + double ub); absl::StatusOr AppendCons(const std::vector& lb, const std::vector& ub); absl::Status PutVarType(VariableIndex j, bool is_integer); absl::Status PutC(const std::vector& c); + absl::Status PutQObj(const std::vector & subk, + const std::vector & subl, + const std::vector & valkl); + absl::Status UpdateQObjEntries(const std::vector & subk, + const std::vector & subl, + const std::vector & valkl); void PutCFix(double cfix); + absl::Status PutQCon(int i, const std::vector& subk, + const std::vector& subl, + const std::vector& cof); + absl::Status PutARow(int i, const std::vector & subj, const std::vector & cof); absl::Status PutAIJList(const std::vector& subi, const std::vector& subj, const std::vector& valij); @@ -173,6 +185,9 @@ class Mosek { static int info_callback(MSKtask_t task, MSKuserhandle_t h, MSKcallbackcodee code, const double* dinf, const int* iinf, const int64_t* liinf); + + absl::StatusOr append_cons(int n, const double* lb, + const double* ub); }; } // namespace operations_research::math_opt diff --git a/ortools/math_opt/solvers/mosek_solver.cc b/ortools/math_opt/solvers/mosek_solver.cc index 60c29f8012d..e3539f0bce4 100644 --- a/ortools/math_opt/solvers/mosek_solver.cc +++ b/ortools/math_opt/solvers/mosek_solver.cc @@ -25,6 +25,7 @@ #include #include #include +#include #include #include @@ -56,6 +57,7 @@ #include "ortools/math_opt/solution.pb.h" #include "ortools/math_opt/solvers/message_callback_data.h" #include "ortools/math_opt/solvers/mosek.pb.h" +#include "ortools/math_opt/sparse_containers.pb.h" #include "ortools/util/solve_interrupter.h" #include "ortools/util/status_macros.h" @@ -68,6 +70,7 @@ constexpr SupportedProblemStructures kMosekSupportedStructures = { .indicator_constraints = SupportType::kSupported, }; + } // namespace std::ostream& mosek::operator<<(std::ostream& s, SolSta solsta) { @@ -198,9 +201,117 @@ absl::Status MosekSolver::ReplaceObjective(const ObjectiveProto& obj) { c[objcof.ids(i)] = objcof.values(i); } RETURN_IF_ERROR(msk.PutC(c)); + + // quadratic terms + if (obj.quadratic_coefficients().row_ids_size() > 0) { + std::vector subk; + std::vector subl; + std::vector val; + + SparseDoubleMatrixToTril(obj.quadratic_coefficients(),subk,subl,val); + RETURN_IF_ERROR(msk.PutQObj(subk,subl,val)); + } + return absl::OkStatus(); } // MosekSolver::ReplaceObjective +void MosekSolver::SparseDoubleMatrixToTril(const SparseDoubleMatrixProto& qdata, + std::vector& subk, + std::vector& subl, + std::vector& val) { + if (qdata.row_ids_size() > 0) { + // NOTE: this specifies the full Q matrix, and we assume that it is + // symmetric and only specifies the lower triangular part. + size_t nqnz = qdata.row_ids_size(); + std::vector> subklv; subklv.reserve(nqnz); + for (auto [kit,lit,cit] = std::make_tuple( + qdata.row_ids().begin(), + qdata.column_ids().begin(), + qdata.coefficients().begin()); + kit != qdata.row_ids().end() && + lit != qdata.column_ids().end() && + cit != qdata.coefficients().end(); + ++kit, ++lit, ++cit) { + if (variable_map.contains(*kit) && + variable_map.contains(*lit)) { + int k = variable_map[*kit]; + int l = variable_map[*lit]; + int v = variable_map[*cit]; + + if (k < l) { + subklv.push_back({l,k,v}); + } + else { + subklv.push_back({k,l,v}); + } + } + } + + std::sort(subklv.begin(),subklv.end()); + + // count + size_t nunique = 0; + { + int prevk = -1,prevl = -1; + for (auto [k,l,v] : subklv) { + if (prevk != k || prevl != l) { + ++nunique; prevk = k; prevl = l; + } + } + } + + subk.clear(); subk.reserve(nunique); + subl.clear(); subl.reserve(nunique); + val.clear(); val.reserve(nunique); + + int prevk = -1,prevl = -1; + for (auto [k,l,v] : subklv) { + if (prevk == k && prevl == l) { + val.back() += v; + } + else { + subk.push_back(k); prevk = k; + subk.push_back(l); prevl = l; + val.push_back(v); + } + } + } +} + +absl::Status MosekSolver::AddConstraint(int64_t id, const QuadraticConstraintProto& cons) { + int coni = msk.NumCon(); + double clb = cons.lower_bound(); + double cub = cons.upper_bound(); + { + auto r = msk.AppendCons(clb, cub); + if (!r.ok()) return r.status(); + } + + size_t nnz = cons.linear_terms().ids_size(); + std::vector subj; subj.reserve(nnz); + std::vector val; val.reserve(nnz); + + for (const auto id : cons.linear_terms().ids()) subj.push_back(variable_map[id]); + for (const auto c : cons.linear_terms().values()) val.push_back(c); + RETURN_IF_ERROR(msk.PutARow(coni,subj, val)); + + // quadratic terms + + if (cons.quadratic_terms().row_ids_size() > 0) { + std::vector subk; + std::vector subl; + std::vector val; + + SparseDoubleMatrixToTril(cons.quadratic_terms(),subk,subl,val); + + RETURN_IF_ERROR(msk.PutQCon(coni,subk,subl,val)); + } + + quadconstr_map[id] = coni; + + return absl::OkStatus(); +} + absl::Status MosekSolver::AddConstraints(const LinearConstraintsProto& cons, const SparseDoubleMatrixProto& adata) { int firstcon = msk.NumCon(); @@ -370,22 +481,46 @@ absl::Status MosekSolver::AddConicConstraints( absl::StatusOr MosekSolver::Update(const ModelUpdateProto& model_update) { for (auto id : model_update.deleted_variable_ids()) { - variable_map.erase(id); - RETURN_IF_ERROR(msk.ClearVariable(variable_map[id])); + if (variable_map.contains(id)) { + int j = variable_map[id]; + variable_map.erase(id); + RETURN_IF_ERROR(msk.ClearVariable(j)); + } } for (auto id : model_update.deleted_linear_constraint_ids()) { - linconstr_map.erase(id); - RETURN_IF_ERROR(msk.ClearConstraint(linconstr_map[id])); + if (linconstr_map.contains(id)) { + int i = linconstr_map[id]; + linconstr_map.erase(id); + RETURN_IF_ERROR(msk.ClearConstraint(i)); + } } for (auto id : model_update.second_order_cone_constraint_updates() .deleted_constraint_ids()) { - coneconstr_map.erase(id); - RETURN_IF_ERROR(msk.ClearConeConstraint(coneconstr_map[id])); + if (coneconstr_map.contains(id)) { + int64_t i = coneconstr_map[id]; + coneconstr_map.erase(id); + RETURN_IF_ERROR(msk.ClearConeConstraint(i)); + } } for (auto id : model_update.indicator_constraint_updates().deleted_constraint_ids()) { - indconstr_map.erase(id); - RETURN_IF_ERROR(msk.ClearDisjunctiveConstraint(indconstr_map[id])); + if (indconstr_map.contains(id)) { + int64_t i = indconstr_map[id]; + indconstr_map.erase(id); + RETURN_IF_ERROR(msk.ClearDisjunctiveConstraint(i)); + } + } + + for (auto id : model_update.quadratic_constraint_updates().deleted_constraint_ids()) { + if (quadconstr_map.contains(id)) { + int i = quadconstr_map[id]; + quadconstr_map.erase(id); + RETURN_IF_ERROR(msk.ClearConstraint(i)); + } + } + + for (auto &[id,con] : model_update.quadratic_constraint_updates().new_constraints()) { + RETURN_IF_ERROR(AddConstraint(id, con)); } RETURN_IF_ERROR(AddVariables(model_update.new_variables())); @@ -469,6 +604,58 @@ absl::Status MosekSolver::UpdateObjective( for (auto id : objupds.linear_coefficients().ids()) subj.push_back(variable_map[id]); + + if (objupds.quadratic_coefficients().column_ids_size() > 0) { + // note: this specifies the full q matrix, and we assume that it is + // symmetric and only specifies the lower triangular part. + auto & qobj = objupds.quadratic_coefficients(); + size_t nqnz = qobj.row_ids_size(); + std::vector> subklv; subklv.reserve(nqnz); + for (auto [kit,lit,cit] = std::make_tuple( + qobj.row_ids().begin(), + qobj.column_ids().begin(), + qobj.coefficients().begin()); + kit != qobj.row_ids().end() && + lit != qobj.column_ids().end() && + cit != qobj.coefficients().end(); + ++kit, ++lit, ++cit) { + int k = variable_map[*kit]; + int l = variable_map[*lit]; + int v = variable_map[*cit]; + + if (k < l) { + subklv.push_back({l,k,v}); + } + else { + subklv.push_back({k,l,v}); + } + } + + std::sort(subklv.begin(),subklv.end()); + + std::vector subk; subk.reserve(nqnz); + std::vector subl; subl.reserve(nqnz); + std::vector val; val.reserve(nqnz); + + int prevk = -1,prevl = -1; + for (auto [k,l,v] : subklv) { + if (prevk == k && prevl == l) { + val.back() += v; + } + else { + subk.push_back(k); prevk = k; + subk.push_back(l); prevl = l; + val.push_back(v); + } + } + + } + + + + + + RETURN_IF_ERROR(msk.UpdateObjectiveSense(objupds.direction_update())); RETURN_IF_ERROR(msk.UpdateObjective(objupds.offset_update(), subj, cof)); diff --git a/ortools/math_opt/solvers/mosek_solver.h b/ortools/math_opt/solvers/mosek_solver.h index 12aeb992169..ff2b2afa1f5 100644 --- a/ortools/math_opt/solvers/mosek_solver.h +++ b/ortools/math_opt/solvers/mosek_solver.h @@ -59,11 +59,13 @@ class MosekSolver : public SolverInterface { Mosek msk; absl::flat_hash_map variable_map; absl::flat_hash_map linconstr_map; + absl::flat_hash_map quadconstr_map; absl::flat_hash_map coneconstr_map; absl::flat_hash_map indconstr_map; absl::Status ReplaceObjective(const ObjectiveProto& obj); absl::Status AddVariables(const VariablesProto& vars); + absl::Status AddConstraint(int64_t id, const QuadraticConstraintProto& cons); absl::Status AddConstraints(const LinearConstraintsProto& cons); absl::Status AddConstraints(const LinearConstraintsProto& cons, const SparseDoubleMatrixProto& lincofupds); @@ -101,86 +103,13 @@ class MosekSolver : public SolverInterface { bool skip_y_zeros, const std::vector& ordered_yx_ids, bool skip_yx_zeros); + void SparseDoubleMatrixToTril(const SparseDoubleMatrixProto & qdata, + std::vector & subi, + std::vector & subj, + std::vector & cof); + MosekSolver(Mosek&& msk); MosekSolver(MosekSolver&) = delete; - -#if 0 - struct SolutionClaims { - bool mosek_returned_primal_feasible_solution = false; - bool mosek_returned_dual_feasible_solution = false; - bool mosek_returned_primal_ray = false; - bool mosek_returned_dual_ray = false; - }; - struct SolutionsAndClaims { - std::vector solutions; - SolutionClaims solution_claims; - }; - - - absl::StatusOr PrimalRayReturned() const; - absl::StatusOr DualRayReturned() const; - - // Will append solutions and rays to result_out. No other modifications are - // made to result_out. Requires that mosek_->getInfo is validated. - absl::StatusOr ExtractSolutionAndRays( - const ModelSolveParametersProto& model_params); - - static absl::StatusOr PrimalFeasibilityStatus( - SolutionClaims solution_claims); - static absl::StatusOr DualFeasibilityStatus( - const MosekInfo& mosek_info, bool is_integer, - SolutionClaims solution_claims); - static absl::StatusOr MakeTermination( - MosekModelStatus mosek_model_status, const MosekInfo& mosek_info, - bool is_integer, bool had_node_limit, bool had_solution_limit, - bool is_maximize, SolutionClaims solution_claims); - - // Returns the current basis if it is available and MathOpt can represent it - // (all kNonBasic values can be made more precise, see b/272767311). - absl::StatusOr> ExtractBasis(); - - template - absl::Status EnsureOneEntryPerVariable(const std::vector& vec) { - if (vec.size() != variable_data_.size()) { - return util::InvalidArgumentErrorBuilder() - << "expected one entry per variable, but model had " - << variable_data_.size() << " variables and found " << vec.size() - << " elements"; - } - return absl::OkStatus(); - } - - template - absl::Status EnsureOneEntryPerLinearConstraint(const std::vector& vec) { - if (vec.size() != lin_con_data_.size()) { - return util::InvalidArgumentErrorBuilder() - << "expected one entry per linear constraint, but model had " - << lin_con_data_.size() << " linear constraints and found " - << vec.size() << " elements"; - } - return absl::OkStatus(); - } - - // Returns a SolveResult for when HiGHS returns the model status - // MosekModelStatus::kModelEmpty. This happens on models that have no - // variables, but may still have (potentially infeasible) linear constraints - // and an objective offset. - // - // Assumes that there are no inverted linear constraint bounds. - static SolveResultProto ResultForMosekModelStatusModelEmpty( - bool is_maximize, double objective_offset, - const absl::flat_hash_map& lin_con_data); - - InvertedBounds ListInvertedBounds(); - - std::unique_ptr mosek_; - - // Key is the mathopt id, value.index is the variable index in HiGHS. - absl::flat_hash_map variable_data_; - - // Key is the mathopt id, value.index is the linear constraint index in HiGHS. - absl::flat_hash_map lin_con_data_; -#endif }; } // namespace operations_research::math_opt From fde4a63b9c55f3493c4cf3adc46cd64eb3e27913 Mon Sep 17 00:00:00 2001 From: ulfw Date: Wed, 30 Oct 2024 12:38:10 +0100 Subject: [PATCH 04/10] working on quadratic support for mosek --- ortools/math_opt/solvers/mosek/mosekwrp.cc | 85 ++++++++++++--- ortools/math_opt/solvers/mosek/mosekwrp.h | 9 +- ortools/math_opt/solvers/mosek_solver.cc | 119 ++++++++++++++------- 3 files changed, 154 insertions(+), 59 deletions(-) diff --git a/ortools/math_opt/solvers/mosek/mosekwrp.cc b/ortools/math_opt/solvers/mosek/mosekwrp.cc index ecb4a8f553d..8cbda1dfff4 100644 --- a/ortools/math_opt/solvers/mosek/mosekwrp.cc +++ b/ortools/math_opt/solvers/mosek/mosekwrp.cc @@ -6,6 +6,7 @@ #include #include #include +#include #include #include "absl/cleanup/cleanup.h" @@ -54,7 +55,10 @@ absl::StatusOr Mosek::AppendVars( if (n > std::numeric_limits::max()) return absl::InvalidArgumentError("arguments lb and ub too large"); - MSK_appendvars(task.get(), (int)n); + if (MSK_RES_OK != MSK_appendvars(task.get(), (int)n)) { + auto [msg,rname,r] = std::move(LastError()); + return absl::InternalError((std::stringstream() << rname << "(" << r << "): " << msg).str()); + } std::vector bk(n); for (ssize_t i = 0; i < n; ++i) { bk[i] = (lb[i] > ub[i] @@ -89,8 +93,10 @@ absl::StatusOr Mosek::append_cons( int n, const double * lb, const double * ub) { int firsti = NumCon(); - if (MSK_RES_OK != MSK_appendcons(task.get(), n)) - return absl::InternalError("Maximum number of internal constraints exceeded"); + if (MSK_RES_OK != MSK_appendcons(task.get(), n)) { + auto [msg,rname,r] = std::move(LastError()); + return absl::InternalError((std::stringstream() << rname << "(" << r << "): " << msg).str()); + } std::vector bk(n); for (ssize_t i = 0; i < n; ++i) { bk[i] = (lb[i] > ub[i] @@ -109,8 +115,9 @@ absl::StatusOr Mosek::append_cons( absl::Status Mosek::PutVarType(VariableIndex j, bool is_integer) { if (MSK_RES_OK != MSK_putvartype(task.get(), j, - is_integer ? MSK_VAR_TYPE_INT : MSK_VAR_TYPE_CONT)) + is_integer ? MSK_VAR_TYPE_INT : MSK_VAR_TYPE_CONT)) { return absl::InvalidArgumentError("Arguments j is invalid"); + } return absl::OkStatus(); } @@ -125,10 +132,12 @@ absl::Status Mosek::PutQObj(const std::vector & subk, int qnnz = (int)subk.size(); if (qnnz != subk.size()) return absl::InvalidArgumentError("Number of quadratic objective nonzeros is too large"); - + if (MSK_RES_OK != - MSK_putqobj(task.get(),qnnz,subk.data(),subl.data(),valkl.data())) - return absl::InternalError("Error setting quadratic objective terms"); + MSK_putqobj(task.get(),qnnz,subk.data(),subl.data(),valkl.data())) { + auto [msg,rname,r] = std::move(LastError()); + return absl::InternalError((std::stringstream() << rname << "(" << r << "): " << msg).str()); + } return absl::OkStatus(); } absl::Status Mosek::UpdateQObjEntries(const std::vector & subk, @@ -146,8 +155,10 @@ absl::Status Mosek::UpdateQObjEntries(const std::vector & subk, ++kit, ++lit, ++cit) r = MSK_putqobjij(task.get(),*kit,*lit,*cit); - if (r != MSK_RES_OK) - return absl::InternalError("Failed to input quadratic objective nonzeros"); + if (r != MSK_RES_OK) { + auto [msg,rname,r] = std::move(LastError()); + return absl::InternalError((std::stringstream() << rname << "(" << r << "): " << msg).str()); + } return absl::OkStatus(); } @@ -168,8 +179,10 @@ absl::Status Mosek::PutQCon(int i, const std::vector& subk, subk.size() != cof.size()) return absl::InvalidArgumentError("Mismatching argument sizes"); if (MSK_RES_OK != - MSK_putqconk(task.get(),i,subk.size(),subk.data(),subl.data(),cof.data())) - return absl::InternalError("Failed to set quadratic nonzeros in constraint"); + MSK_putqconk(task.get(),i,subk.size(),subk.data(),subl.data(),cof.data())) { + auto [msg,rname,r] = std::move(LastError()); + return absl::InternalError((std::stringstream() << rname << "(" << r << "): " << msg).str()); + } return absl::OkStatus(); } @@ -182,10 +195,13 @@ absl::Status Mosek::PutARow(int i, const std::vector& subj, return absl::InvalidArgumentError("Number of nonzeros added too large"); int nnz = (int)subj.size(); - if (MSK_RES_OK != MSK_putarow(task.get(),i,nnz,subj.data(),cof.data())) - return absl::InternalError("Failed to add row nonzeros"); + if (MSK_RES_OK != MSK_putarow(task.get(),i,nnz,subj.data(),cof.data())) { + auto [msg,rname,r] = std::move(LastError()); + return absl::InternalError((std::stringstream() << rname << "(" << r << "): " << msg).str()); + } return absl::OkStatus(); } + absl::Status Mosek::PutAIJList(const std::vector& subi, const std::vector& subj, const std::vector& valij) { @@ -198,8 +214,10 @@ absl::Status Mosek::PutAIJList(const std::vector& subi, "Arguments subi, subj, valij are too long"); if (MSK_RES_OK != MSK_putaijlist(task.get(), (int)n, subi.data(), subj.data(), - valij.data())) - return absl::InvalidArgumentError("Invalid index argument subi or subj"); + valij.data())) { + auto [msg,rname,r] = std::move(LastError()); + return absl::InternalError((std::stringstream() << rname << "(" << r << "): " << msg).str()); + } return absl::OkStatus(); } @@ -462,7 +480,10 @@ void Mosek::WriteData(const char* filename) const { absl::StatusOr Mosek::Optimize() { MSKrescodee trm; MSKrescodee r = MSK_optimizetrm(task.get(), &trm); - if (MSK_RES_OK != r) return absl::InternalError("Optimization failed"); + if (MSK_RES_OK != r) { + std::cerr << "MOSEK ERROR: " << r << std::endl; + return absl::InternalError("Optimization failedcode"); + } return trm; } @@ -501,7 +522,10 @@ absl::StatusOr Mosek::Optimize(MessageCallback msg_cb, message_callback); MSKrescodee r = MSK_optimizetrm(task.get(), &trm); - if (MSK_RES_OK != r) return absl::InternalError("Optimization failed"); + if (MSK_RES_OK != r) { + std::cerr << "MOSEK ERROR: " << r << std::endl; + return absl::InternalError("Optimization failed"); + } return trm; } @@ -522,6 +546,8 @@ std::tuple Mosek::LastError() const { } } + + double Mosek::GetPrimalObj(MSKsoltypee whichsol) const { if (!SolutionDef(whichsol)) return 0.0; double val; @@ -583,6 +609,22 @@ void Mosek::GetSKC(MSKsoltypee whichsol, std::vector& skc) const { MSK_getskc(task.get(), whichsol, skc.data()); } +int Mosek::GetIntInfoItem(MSKiinfiteme item) { + int val = 0; + MSK_getintinf(task.get(),item,&val); + return val; +} +int64_t Mosek::GetLongInfoItem(MSKliinfiteme item) { + int64_t val = 0; + MSK_getlintinf(task.get(),item,&val); + return val; +} +double Mosek::GetDoubleInfoItem(MSKdinfiteme item) { + double val = 0; + MSK_getdouinf(task.get(),item,&val); + return val; +} + void Mosek::PutParam(MSKiparame par, int value) { MSK_putintparam(task.get(), par, value); } @@ -600,6 +642,13 @@ int Mosek::NumCon() const { MSK_getnumcon(task.get(), &n); return n; } + +void Mosek::GetC(std::vector & c, double & cfix) { + c.resize(NumVar()); + MSK_getc(task.get(),c.data()); + MSK_getcfix(task.get(),&cfix); +} + bool Mosek::IsMaximize() const { MSKobjsensee sense; MSK_getobjsense(task.get(), &sense); @@ -617,4 +666,6 @@ int Mosek::GetParam(MSKiparame ipar) const { } void Mosek::delete_msk_task_func(MSKtask_t t) { MSK_deletetask(&t); } + + } // namespace operations_research::math_opt diff --git a/ortools/math_opt/solvers/mosek/mosekwrp.h b/ortools/math_opt/solvers/mosek/mosekwrp.h index dbde1eb4bb8..70fb5f7d644 100644 --- a/ortools/math_opt/solvers/mosek/mosekwrp.h +++ b/ortools/math_opt/solvers/mosek/mosekwrp.h @@ -141,6 +141,8 @@ class Mosek { double GetParam(MSKdparame dpar) const; int GetParam(MSKiparame ipar) const; + void GetC(std::vector & c, double & cfix); + bool SolutionDef(MSKsoltypee which) const { MSKbooleant soldef; MSK_solutiondef(task.get(), which, &soldef); @@ -157,6 +159,11 @@ class Mosek { return (mosek::SolSta)solsta; } + + int GetIntInfoItem(MSKiinfiteme item); + int64_t GetLongInfoItem(MSKliinfiteme item); + double GetDoubleInfoItem(MSKdinfiteme item); + std::tuple LastError() const; double GetPrimalObj(MSKsoltypee whichsol) const; @@ -185,7 +192,7 @@ class Mosek { static int info_callback(MSKtask_t task, MSKuserhandle_t h, MSKcallbackcodee code, const double* dinf, const int* iinf, const int64_t* liinf); - + absl::StatusOr append_cons(int n, const double* lb, const double* ub); }; diff --git a/ortools/math_opt/solvers/mosek_solver.cc b/ortools/math_opt/solvers/mosek_solver.cc index e3539f0bce4..9cd74b77fea 100644 --- a/ortools/math_opt/solvers/mosek_solver.cc +++ b/ortools/math_opt/solvers/mosek_solver.cc @@ -1,4 +1,4 @@ -// Copyright 2010-2024 Google LLskip_xx_zeros C +// Copyright 2010-2024 Mosek // 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 @@ -19,38 +19,37 @@ #include #include #include -#include +#include +//#include +//#include #include #include #include -#include +//#include #include #include #include #include -// #include "absl/algorithm/container.h" -// #include "absl/cleanup/cleanup.h" -// #include "absl/container/flat_hash_map.h" #include "absl/log/check.h" -// #include "absl/memory/memory.h" +#include "absl/time/clock.h" +#include "absl/time/time.h" #include "absl/cleanup/cleanup.h" #include "absl/status/status.h" #include "absl/status/statusor.h" #include "absl/strings/string_view.h" -// #include "absl/time/clock.h" -// #include "absl/time/time.h" #include "absl/types/span.h" +#include "mosek.h" #include "ortools/base/protoutil.h" #include "ortools/base/status_builder.h" #include "ortools/base/status_macros.h" #include "ortools/math_opt/callback.pb.h" -#include "ortools/math_opt/core/empty_bounds.h" -#include "ortools/math_opt/core/inverted_bounds.h" +//#include "ortools/math_opt/core/empty_bounds.h" +//#include "ortools/math_opt/core/inverted_bounds.h" #include "ortools/math_opt/core/math_opt_proto_utils.h" #include "ortools/math_opt/core/solver_interface.h" -#include "ortools/math_opt/core/sorted.h" -#include "ortools/math_opt/core/sparse_vector_view.h" +//#include "ortools/math_opt/core/sorted.h" +//#include "ortools/math_opt/core/sparse_vector_view.h" #include "ortools/math_opt/infeasible_subsystem.pb.h" #include "ortools/math_opt/parameters.pb.h" #include "ortools/math_opt/result.pb.h" @@ -66,6 +65,8 @@ namespace { constexpr SupportedProblemStructures kMosekSupportedStructures = { .integer_variables = SupportType::kSupported, + .quadratic_objectives = SupportType::kSupported, + .quadratic_constraints = SupportType::kSupported, .second_order_cone_constraints = SupportType::kSupported, .indicator_constraints = SupportType::kSupported, }; @@ -209,6 +210,11 @@ absl::Status MosekSolver::ReplaceObjective(const ObjectiveProto& obj) { std::vector val; SparseDoubleMatrixToTril(obj.quadratic_coefficients(),subk,subl,val); + //std::cerr << "subk,subl,val: " << subk.size() << "," << subl.size() << "," << val.size() << std::endl; + //for (auto [ik,il,iv] = std::make_tuple(subk.begin(),subl.begin(),val.begin()); + // ik != subk.end(); + // ++ik,++il,++iv) + // std::cerr << " -- [" << *ik << "," << *il << "]: " << *iv << std::endl; RETURN_IF_ERROR(msk.PutQObj(subk,subl,val)); } @@ -224,19 +230,19 @@ void MosekSolver::SparseDoubleMatrixToTril(const SparseDoubleMatrixProto& qdata, // symmetric and only specifies the lower triangular part. size_t nqnz = qdata.row_ids_size(); std::vector> subklv; subklv.reserve(nqnz); - for (auto [kit,lit,cit] = std::make_tuple( - qdata.row_ids().begin(), - qdata.column_ids().begin(), - qdata.coefficients().begin()); - kit != qdata.row_ids().end() && - lit != qdata.column_ids().end() && + for (auto [kit, lit, cit] = std::make_tuple(qdata.row_ids().begin(), + qdata.column_ids().begin(), + qdata.coefficients().begin()); + kit != qdata.row_ids().end() && lit != qdata.column_ids().end() && cit != qdata.coefficients().end(); ++kit, ++lit, ++cit) { + if (variable_map.contains(*kit) && variable_map.contains(*lit)) { int k = variable_map[*kit]; int l = variable_map[*lit]; - int v = variable_map[*cit]; + double v = *cit; + if (k < l) { subklv.push_back({l,k,v}); @@ -271,7 +277,7 @@ void MosekSolver::SparseDoubleMatrixToTril(const SparseDoubleMatrixProto& qdata, } else { subk.push_back(k); prevk = k; - subk.push_back(l); prevl = l; + subl.push_back(l); prevl = l; val.push_back(v); } } @@ -690,14 +696,14 @@ absl::StatusOr> MosekSolver::New( if (!model.auxiliary_objectives().empty()) return util::InvalidArgumentErrorBuilder() << "Mosek does not support multi-objective models"; - if (!model.objective().quadratic_coefficients().row_ids().empty()) { - return util::InvalidArgumentErrorBuilder() - << "Mosek does not support models with quadratic objectives"; - } - if (!model.quadratic_constraints().empty()) { - return util::InvalidArgumentErrorBuilder() - << "Mosek does not support models with quadratic constraints"; - } + //if (!model.objective().quadratic_coefficients().row_ids().empty()) { + // return util::InvalidArgumentErrorBuilder() + // << "Mosek does not support models with quadratic objectives"; + //} + //if (!model.quadratic_constraints().empty()) { + // return util::InvalidArgumentErrorBuilder() + // << "Mosek does not support models with quadratic constraints"; + // } if (!model.sos1_constraints().empty() || !model.sos2_constraints().empty()) { return util::InvalidArgumentErrorBuilder() << "Mosek does not support models with SOS constraints"; @@ -711,6 +717,8 @@ absl::StatusOr> MosekSolver::New( RETURN_IF_ERROR(mskslv->ReplaceObjective(model.objective())); RETURN_IF_ERROR(mskslv->AddConstraints(model.linear_constraints(), model.linear_constraint_matrix())); + for (auto &[k,v] : model.quadratic_constraints()) + RETURN_IF_ERROR(mskslv->AddConstraint(k,v)); RETURN_IF_ERROR( mskslv->AddIndicatorConstraints(model.indicator_constraints())); @@ -990,6 +998,9 @@ absl::StatusOr MosekSolver::Solve( // - EmphasisProto cuts // - EmphasisProto heuristics // - EmphasisProto scaling + + + auto solve_start = absl::Now(); double dpar_optimizer_max_time = msk.GetParam(MSK_DPAR_OPTIMIZER_MAX_TIME); int ipar_intpnt_max_iterations = msk.GetParam(MSK_IPAR_INTPNT_MAX_ITERATIONS); @@ -1138,10 +1149,11 @@ absl::StatusOr MosekSolver::Solve( ordered_yx_ids.push_back(id); } - MSKrescodee trm; + MSKrescodee trm = MSK_RES_OK; { BufferedMessageCallback bmsg_cb(message_cb); // TODO: Use model_parameters + msk.WriteData("test.opf"); auto r = msk.Optimize( [&](const std::string& msg) { bmsg_cb.OnMessage(msg); }, [&](MSKcallbackcodee code, absl::Span dinf, @@ -1246,8 +1258,6 @@ absl::StatusOr MosekSolver::Solve( msk.IsMaximize(), TerminationReasonProto::TERMINATION_REASON_NO_SOLUTION_FOUND, msg); } else { - // std::cout << "MosekSolver::Solve() solution is defined " << whichsol << - // std::endl; prosta = msk.GetProSta(whichsol); solsta = msk.GetSolSta(whichsol); @@ -1256,9 +1266,15 @@ absl::StatusOr MosekSolver::Solve( if (solsta == mosek::SolSta::OPTIMAL || solsta == mosek::SolSta::INTEGER_OPTIMAL) { - // std::cout << "MosekSolver::Solve() trmp = Optimal! " << std::endl; + trmp = OptimalTerminationProto(msk.GetPrimalObj(whichsol), - msk.GetDualObj(whichsol), ""); + msk.GetDualObj(whichsol), + ""); + // Hack: + double dx = msk.IsMaximize() ? -1e-9 : 1e-9; + trmp.mutable_objective_bounds()->set_primal_bound(trmp.objective_bounds().primal_bound()-dx); + trmp.mutable_objective_bounds()->set_dual_bound(trmp.objective_bounds().dual_bound()+dx); + } else if (solsta == mosek::SolSta::PRIM_INFEAS_CER) trmp = InfeasibleTerminationProto( msk.IsMaximize(), @@ -1325,14 +1341,24 @@ absl::StatusOr MosekSolver::Solve( case mosek::SolSta::INTEGER_OPTIMAL: case mosek::SolSta::PRIM_FEAS: case mosek::SolSta::DUAL_FEAS: - case mosek::SolSta::PRIM_AND_DUAL_FEAS: { - auto r = Solution(whichsol, ordered_xc_ids, ordered_xx_ids, - skip_xx_zeros, ordered_y_ids, skip_y_zeros, - ordered_yx_ids, skip_yx_zeros); - if (r.ok()) { - *result.add_solutions() = std::move(*r); + case mosek::SolSta::PRIM_AND_DUAL_FEAS: + { + auto r = Solution(whichsol, ordered_xc_ids, ordered_xx_ids, + skip_xx_zeros, ordered_y_ids, skip_y_zeros, + ordered_yx_ids, skip_yx_zeros); + if (r.ok()) { + *result.add_solutions() = std::move(*r); + } + //if (msk.IsMaximize()) { + // result.mutable_termination()->mutable_objective_bounds()->set_primal_bound(- std::numeric_limits::infinity()); + // result.mutable_termination()->mutable_objective_bounds()->set_dual_bound(std::numeric_limits::infinity()); + //} + //else { + // result.mutable_termination()->mutable_objective_bounds()->set_primal_bound(std::numeric_limits::infinity()); + // result.mutable_termination()->mutable_objective_bounds()->set_dual_bound(- std::numeric_limits::infinity()); + //} } - } break; + break; case mosek::SolSta::DUAL_INFEAS_CER: { auto r = PrimalRay(whichsol, ordered_xx_ids, skip_xx_zeros); if (r.ok()) { @@ -1352,6 +1378,17 @@ absl::StatusOr MosekSolver::Solve( break; } } + + SolveStatsProto * stats = result.mutable_solve_stats(); + stats->set_simplex_iterations(msk.GetIntInfoItem(MSK_IINF_SIM_PRIMAL_ITER)+msk.GetIntInfoItem(MSK_IINF_SIM_DUAL_ITER)); + stats->set_barrier_iterations(msk.GetIntInfoItem(MSK_IINF_INTPNT_ITER)); + stats->set_node_count(msk.GetIntInfoItem(MSK_IINF_MIO_NUM_SOLVED_NODES)); + + auto solve_time_proto = util_time::EncodeGoogleApiProto(absl::Now() - solve_start); + if (solve_time_proto.ok()) { + *stats->mutable_solve_time() = *solve_time_proto; + } + return result; } From a470460f1cc59602afb9843c8c8717e53ca5b02c Mon Sep 17 00:00:00 2001 From: ulfw Date: Wed, 30 Oct 2024 12:55:25 +0100 Subject: [PATCH 05/10] fix mosek primal/dual obj mixup --- ortools/math_opt/solvers/mosek_solver.cc | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/ortools/math_opt/solvers/mosek_solver.cc b/ortools/math_opt/solvers/mosek_solver.cc index 9cd74b77fea..ae180eb565b 100644 --- a/ortools/math_opt/solvers/mosek_solver.cc +++ b/ortools/math_opt/solvers/mosek_solver.cc @@ -772,7 +772,7 @@ absl::StatusOr MosekSolver::DualSolution( case mosek::SolSta::OPTIMAL: case mosek::SolSta::PRIM_AND_DUAL_FEAS: case mosek::SolSta::DUAL_FEAS: { - sol.set_objective_value(msk.GetPrimalObj(whichsol)); + sol.set_objective_value(msk.GetDualObj(whichsol)); sol.set_feasibility_status(SolutionStatusProto::SOLUTION_STATUS_FEASIBLE); std::vector keys; keys.reserve(std::max(variable_map.size(), linconstr_map.size())); @@ -1271,9 +1271,8 @@ absl::StatusOr MosekSolver::Solve( msk.GetDualObj(whichsol), ""); // Hack: - double dx = msk.IsMaximize() ? -1e-9 : 1e-9; - trmp.mutable_objective_bounds()->set_primal_bound(trmp.objective_bounds().primal_bound()-dx); - trmp.mutable_objective_bounds()->set_dual_bound(trmp.objective_bounds().dual_bound()+dx); + trmp.mutable_objective_bounds()->set_primal_bound(trmp.objective_bounds().primal_bound()); + trmp.mutable_objective_bounds()->set_dual_bound(trmp.objective_bounds().dual_bound()); } else if (solsta == mosek::SolSta::PRIM_INFEAS_CER) trmp = InfeasibleTerminationProto( @@ -1345,8 +1344,8 @@ absl::StatusOr MosekSolver::Solve( { auto r = Solution(whichsol, ordered_xc_ids, ordered_xx_ids, skip_xx_zeros, ordered_y_ids, skip_y_zeros, - ordered_yx_ids, skip_yx_zeros); - if (r.ok()) { + ordered_yx_ids, skip_yx_zeros); + if (r.ok()) { *result.add_solutions() = std::move(*r); } //if (msk.IsMaximize()) { From e978115630abb6a42d2baca7044431a3049dc6a0 Mon Sep 17 00:00:00 2001 From: Ulfw Date: Wed, 30 Oct 2024 12:57:49 +0100 Subject: [PATCH 06/10] add example --- examples/cpp/portfolio_1_basic.cc | 100 ++++++++++++++++++++++++++++++ 1 file changed, 100 insertions(+) create mode 100644 examples/cpp/portfolio_1_basic.cc diff --git a/examples/cpp/portfolio_1_basic.cc b/examples/cpp/portfolio_1_basic.cc new file mode 100644 index 00000000000..0ed83776f28 --- /dev/null +++ b/examples/cpp/portfolio_1_basic.cc @@ -0,0 +1,100 @@ +#include +#include +#include +#include +#include + +#include "absl/log/check.h" +#include "absl/memory/memory.h" +#include "absl/status/statusor.h" +#include "ortools/base/init_google.h" +#include "ortools/math_opt/cpp/math_opt.h" + +const double inf = std::numeric_limits::infinity(); +namespace math_opt = ::operations_research::math_opt; +using math_opt::LinearExpression; + + +/// Build and solve a basic Markowit portfolio problem. +/// +/// # Arguments +/// - mu The vector of length n of expected returns on the assets +/// - GT A vector defining a m x n matrix in row-major format. It is the +/// facored co-variance matrix, i.e. the covariance matrix is Q = G*G^T +/// - x0 A vector of length n of initial investments +/// - w Initial wealth not invested yet. +/// - gamma The risk bound as a bound on the variance of the return of portfolio, +/// gamma > sqrt( xGG^Tx ). +absl::StatusOr>> basic_markowitz( + const std::vector & mu, + const std::vector & GT, + const std::vector & x0, + double w, + double gamma) { + + math_opt::Model model("portfolio_1_basic"); + size_t n = mu.size(); + size_t m = GT.size()/n; + std::vector x; x.reserve(n); + for (int i = 0; i < n; ++i) + x.push_back(model.AddContinuousVariable(0.0,inf,(std::ostringstream() << "x" << i).str())); + + model.Maximize(LinearExpression::InnerProduct(x,mu)); + + model.AddLinearConstraint(LinearExpression::Sum(x) == 0.0, "Budget"); + + std::vector linear_to_norm; + for (int i = 0; i < m; ++i) { + linear_to_norm.push_back(LinearExpression::InnerProduct(absl::Span(GT.data()+m*i,m), x)); + } + model.AddSecondOrderConeConstraint(linear_to_norm, gamma, "risk"); + + // Set parameters, e.g. turn on logging. + math_opt::SolveArguments args; + args.parameters.enable_output = true; + + // Solve and ensure an optimal solution was found with no errors. + const absl::StatusOr result = + math_opt::Solve(model, math_opt::SolverType::kMosek, args); + if (! result.ok()) + return result.status(); + + double pobj = result->objective_value(); + std::vector res(n); + for (int i = 0; i < n; ++i) + res[i] = result->variable_values().at(x[i]); + + return std::make_tuple(pobj,std::move(res)); +} + +int main(int argc, char ** argv) { + + const int n = 8; + const double w = 59.0; + std::vector mu = {0.07197349, 0.15518171, 0.17535435, 0.0898094 , 0.42895777, 0.39291844, 0.32170722, 0.18378628}; + std::vector x0 = {8.0, 5.0, 3.0, 5.0, 2.0, 9.0, 3.0, 6.0}; + double gamma = 36; + std::vector GT = { + 0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638, + 0. , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506, + 0. , 0. , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914, + 0. , 0. , 0. , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742, + 0. , 0. , 0. , 0. , 0.36096, 0.12574, 0.10157, 0.0571 , + 0. , 0. , 0. , 0. , 0. , 0.21552, 0.05663, 0.06187, + 0. , 0. , 0. , 0. , 0. , 0. , 0.22514, 0.03327, + 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2202 }; + + + auto res = basic_markowitz(mu,GT,x0,w,gamma); + + if (! res.ok()) { + std::cerr << "Failed to solve problem" << std::endl; + } + else { + auto &[pobj,xx] = *res; + std::cout << "Primal Objective value: " << pobj << std::endl; + std::cout << "Solution values:" << std::endl; + for (int i = 0; i < xx.size(); ++i) + std::cout << " x["<< i << "] = " << xx[i] << std::endl; + } +} From 1b3535c6b2fb38ffdf9811787c1078ee1689e27c Mon Sep 17 00:00:00 2001 From: ulfw Date: Wed, 30 Oct 2024 13:48:17 +0100 Subject: [PATCH 07/10] fix: mosek issues inputting conic constraints --- examples/cpp/portfolio_1_basic.cc | 13 +++---- ortools/math_opt/solvers/mosek/mosekwrp.cc | 8 ++--- ortools/math_opt/solvers/mosek_solver.cc | 41 +++++++++------------- 3 files changed, 28 insertions(+), 34 deletions(-) diff --git a/examples/cpp/portfolio_1_basic.cc b/examples/cpp/portfolio_1_basic.cc index 0ed83776f28..3b3df9919c3 100644 --- a/examples/cpp/portfolio_1_basic.cc +++ b/examples/cpp/portfolio_1_basic.cc @@ -30,7 +30,8 @@ absl::StatusOr>> basic_markowitz( const std::vector & GT, const std::vector & x0, double w, - double gamma) { + double gamma, + math_opt::SolverType st) { math_opt::Model model("portfolio_1_basic"); size_t n = mu.size(); @@ -41,7 +42,8 @@ absl::StatusOr>> basic_markowitz( model.Maximize(LinearExpression::InnerProduct(x,mu)); - model.AddLinearConstraint(LinearExpression::Sum(x) == 0.0, "Budget"); + double totalwealth = w; for (double v : x0) w += v; + model.AddLinearConstraint(LinearExpression::Sum(x) == totalwealth, "Budget"); std::vector linear_to_norm; for (int i = 0; i < m; ++i) { @@ -55,7 +57,7 @@ absl::StatusOr>> basic_markowitz( // Solve and ensure an optimal solution was found with no errors. const absl::StatusOr result = - math_opt::Solve(model, math_opt::SolverType::kMosek, args); + math_opt::Solve(model, st, args); if (! result.ok()) return result.status(); @@ -68,7 +70,6 @@ absl::StatusOr>> basic_markowitz( } int main(int argc, char ** argv) { - const int n = 8; const double w = 59.0; std::vector mu = {0.07197349, 0.15518171, 0.17535435, 0.0898094 , 0.42895777, 0.39291844, 0.32170722, 0.18378628}; @@ -85,10 +86,10 @@ int main(int argc, char ** argv) { 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2202 }; - auto res = basic_markowitz(mu,GT,x0,w,gamma); + auto res = basic_markowitz(mu,GT,x0,w,gamma,math_opt::SolverType::kMosek); if (! res.ok()) { - std::cerr << "Failed to solve problem" << std::endl; + std::cerr << "Failed to solve problem: " << res.status() << std::endl; } else { auto &[pobj,xx] = *res; diff --git a/ortools/math_opt/solvers/mosek/mosekwrp.cc b/ortools/math_opt/solvers/mosek/mosekwrp.cc index 8cbda1dfff4..85e51a2f48e 100644 --- a/ortools/math_opt/solvers/mosek/mosekwrp.cc +++ b/ortools/math_opt/solvers/mosek/mosekwrp.cc @@ -4,10 +4,8 @@ #include #include -#include #include #include -#include #include "absl/cleanup/cleanup.h" #include "absl/status/status.h" @@ -294,8 +292,10 @@ absl::StatusOr Mosek::AppendConeConstraint( if (nnz != cof.size() || nnz != subj.size()) return absl::InvalidArgumentError( "Mismatching argument lengths of subj and cof"); - if (n != b.size()) - return absl::InvalidArgumentError("Mismatching argument lengths b and ptr"); + if (n != b.size()) { + std::cout << " -- b.size() = " << b.size() << ", sizes.size() = " << sizes.size() << std::endl; + return absl::InvalidArgumentError("Mismatching argument lengths b and sizes"); + } switch (ct) { case ConeType::SecondOrderCone: diff --git a/ortools/math_opt/solvers/mosek_solver.cc b/ortools/math_opt/solvers/mosek_solver.cc index ae180eb565b..947953ab7a3 100644 --- a/ortools/math_opt/solvers/mosek_solver.cc +++ b/ortools/math_opt/solvers/mosek_solver.cc @@ -446,6 +446,11 @@ absl::Status MosekSolver::AddConicConstraints( sizes.reserve(cons.size()); for (const auto& [idx, con] : cons) { + subj.clear(); + cof.clear(); + sizes.clear(); + b.clear(); + auto& expr0 = con.upper_bound(); int64_t totalnnz = expr0.ids_size(); for (const auto& lexp : con.arguments_to_norm()) { @@ -454,14 +459,17 @@ absl::Status MosekSolver::AddConicConstraints( subj.reserve(totalnnz); cof.reserve(totalnnz); + b.push_back(expr0.offset()); + for (const auto& expri : con.arguments_to_norm()) { + b.push_back(expri.offset()); + } - for (const auto& id : expr0.ids()) { + sizes.push_back(expr0.ids_size()); + for (const auto& id : expr0.ids()) subj.push_back(variable_map[id]); - } - for (auto c : expr0.coefficients()) { + for (auto c : expr0.coefficients()) cof.push_back(c); - } for (const auto& expri : con.arguments_to_norm()) { sizes.push_back(expri.ids_size()); @@ -471,7 +479,6 @@ absl::Status MosekSolver::AddConicConstraints( for (auto c : expri.coefficients()) { cof.push_back(c); } - b.push_back(expri.offset()); } auto acci = msk.AppendConeConstraint(Mosek::ConeType::SecondOrderCone, @@ -696,14 +703,6 @@ absl::StatusOr> MosekSolver::New( if (!model.auxiliary_objectives().empty()) return util::InvalidArgumentErrorBuilder() << "Mosek does not support multi-objective models"; - //if (!model.objective().quadratic_coefficients().row_ids().empty()) { - // return util::InvalidArgumentErrorBuilder() - // << "Mosek does not support models with quadratic objectives"; - //} - //if (!model.quadratic_constraints().empty()) { - // return util::InvalidArgumentErrorBuilder() - // << "Mosek does not support models with quadratic constraints"; - // } if (!model.sos1_constraints().empty() || !model.sos2_constraints().empty()) { return util::InvalidArgumentErrorBuilder() << "Mosek does not support models with SOS constraints"; @@ -717,6 +716,7 @@ absl::StatusOr> MosekSolver::New( RETURN_IF_ERROR(mskslv->ReplaceObjective(model.objective())); RETURN_IF_ERROR(mskslv->AddConstraints(model.linear_constraints(), model.linear_constraint_matrix())); + RETURN_IF_ERROR(mskslv->AddConicConstraints(model.second_order_cone_constraints())); for (auto &[k,v] : model.quadratic_constraints()) RETURN_IF_ERROR(mskslv->AddConstraint(k,v)); RETURN_IF_ERROR( @@ -724,6 +724,8 @@ absl::StatusOr> MosekSolver::New( std::unique_ptr res(std::move(mskslv)); + //std::cout << " ---------- MosekSolver::New" << std::endl; + return res; } @@ -999,7 +1001,6 @@ absl::StatusOr MosekSolver::Solve( // - EmphasisProto heuristics // - EmphasisProto scaling - auto solve_start = absl::Now(); double dpar_optimizer_max_time = msk.GetParam(MSK_DPAR_OPTIMIZER_MAX_TIME); @@ -1153,7 +1154,7 @@ absl::StatusOr MosekSolver::Solve( { BufferedMessageCallback bmsg_cb(message_cb); // TODO: Use model_parameters - msk.WriteData("test.opf"); + //msk.WriteData("test.ptf"); auto r = msk.Optimize( [&](const std::string& msg) { bmsg_cb.OnMessage(msg); }, [&](MSKcallbackcodee code, absl::Span dinf, @@ -1348,14 +1349,6 @@ absl::StatusOr MosekSolver::Solve( if (r.ok()) { *result.add_solutions() = std::move(*r); } - //if (msk.IsMaximize()) { - // result.mutable_termination()->mutable_objective_bounds()->set_primal_bound(- std::numeric_limits::infinity()); - // result.mutable_termination()->mutable_objective_bounds()->set_dual_bound(std::numeric_limits::infinity()); - //} - //else { - // result.mutable_termination()->mutable_objective_bounds()->set_primal_bound(std::numeric_limits::infinity()); - // result.mutable_termination()->mutable_objective_bounds()->set_dual_bound(- std::numeric_limits::infinity()); - //} } break; case mosek::SolSta::DUAL_INFEAS_CER: { @@ -1401,4 +1394,4 @@ MosekSolver::ComputeInfeasibleSubsystem(const SolveParametersProto&, MATH_OPT_REGISTER_SOLVER(SOLVER_TYPE_MOSEK, MosekSolver::New); -} // namespace operations_research::math_opt +} // namespace operations_research::math_opt From 08904f628d4dad5eff4ea11573fdab966fe43922 Mon Sep 17 00:00:00 2001 From: ulfw Date: Wed, 30 Oct 2024 13:53:10 +0100 Subject: [PATCH 08/10] apply clang format to mosek files --- ortools/math_opt/solvers/mosek_solver.cc | 229 +++++++++++------------ ortools/math_opt/solvers/mosek_solver.h | 7 +- 2 files changed, 116 insertions(+), 120 deletions(-) diff --git a/ortools/math_opt/solvers/mosek_solver.cc b/ortools/math_opt/solvers/mosek_solver.cc index 947953ab7a3..589d1df56b2 100644 --- a/ortools/math_opt/solvers/mosek_solver.cc +++ b/ortools/math_opt/solvers/mosek_solver.cc @@ -20,36 +20,29 @@ #include #include #include -//#include -//#include #include #include #include -//#include #include #include #include #include -#include "absl/log/check.h" -#include "absl/time/clock.h" -#include "absl/time/time.h" #include "absl/cleanup/cleanup.h" +#include "absl/log/check.h" #include "absl/status/status.h" #include "absl/status/statusor.h" #include "absl/strings/string_view.h" +#include "absl/time/clock.h" +#include "absl/time/time.h" #include "absl/types/span.h" #include "mosek.h" #include "ortools/base/protoutil.h" #include "ortools/base/status_builder.h" #include "ortools/base/status_macros.h" #include "ortools/math_opt/callback.pb.h" -//#include "ortools/math_opt/core/empty_bounds.h" -//#include "ortools/math_opt/core/inverted_bounds.h" #include "ortools/math_opt/core/math_opt_proto_utils.h" #include "ortools/math_opt/core/solver_interface.h" -//#include "ortools/math_opt/core/sorted.h" -//#include "ortools/math_opt/core/sparse_vector_view.h" #include "ortools/math_opt/infeasible_subsystem.pb.h" #include "ortools/math_opt/parameters.pb.h" #include "ortools/math_opt/result.pb.h" @@ -71,7 +64,6 @@ constexpr SupportedProblemStructures kMosekSupportedStructures = { .indicator_constraints = SupportType::kSupported, }; - } // namespace std::ostream& mosek::operator<<(std::ostream& s, SolSta solsta) { @@ -203,88 +195,94 @@ absl::Status MosekSolver::ReplaceObjective(const ObjectiveProto& obj) { } RETURN_IF_ERROR(msk.PutC(c)); - // quadratic terms + // quadratic terms if (obj.quadratic_coefficients().row_ids_size() > 0) { std::vector subk; std::vector subl; std::vector val; - SparseDoubleMatrixToTril(obj.quadratic_coefficients(),subk,subl,val); - //std::cerr << "subk,subl,val: " << subk.size() << "," << subl.size() << "," << val.size() << std::endl; - //for (auto [ik,il,iv] = std::make_tuple(subk.begin(),subl.begin(),val.begin()); - // ik != subk.end(); - // ++ik,++il,++iv) - // std::cerr << " -- [" << *ik << "," << *il << "]: " << *iv << std::endl; - RETURN_IF_ERROR(msk.PutQObj(subk,subl,val)); + SparseDoubleMatrixToTril(obj.quadratic_coefficients(), subk, subl, val); + // std::cerr << "subk,subl,val: " << subk.size() << "," << subl.size() << + // "," << val.size() << std::endl; for (auto [ik,il,iv] = + // std::make_tuple(subk.begin(),subl.begin(),val.begin()); + // ik != subk.end(); + // ++ik,++il,++iv) + // std::cerr << " -- [" << *ik << "," << *il << "]: " << *iv << + // std::endl; + RETURN_IF_ERROR(msk.PutQObj(subk, subl, val)); } return absl::OkStatus(); } // MosekSolver::ReplaceObjective void MosekSolver::SparseDoubleMatrixToTril(const SparseDoubleMatrixProto& qdata, - std::vector& subk, - std::vector& subl, - std::vector& val) { + std::vector& subk, + std::vector& subl, + std::vector& val) { if (qdata.row_ids_size() > 0) { // NOTE: this specifies the full Q matrix, and we assume that it is // symmetric and only specifies the lower triangular part. size_t nqnz = qdata.row_ids_size(); - std::vector> subklv; subklv.reserve(nqnz); + std::vector> subklv; + subklv.reserve(nqnz); for (auto [kit, lit, cit] = std::make_tuple(qdata.row_ids().begin(), qdata.column_ids().begin(), qdata.coefficients().begin()); kit != qdata.row_ids().end() && lit != qdata.column_ids().end() && cit != qdata.coefficients().end(); ++kit, ++lit, ++cit) { - - if (variable_map.contains(*kit) && - variable_map.contains(*lit)) { + if (variable_map.contains(*kit) && variable_map.contains(*lit)) { int k = variable_map[*kit]; int l = variable_map[*lit]; double v = *cit; - if (k < l) { - subklv.push_back({l,k,v}); - } - else { - subklv.push_back({k,l,v}); + subklv.push_back({l, k, v}); + } else { + subklv.push_back({k, l, v}); } } } - std::sort(subklv.begin(),subklv.end()); + std::sort(subklv.begin(), subklv.end()); - // count - size_t nunique = 0; + // count + size_t nunique = 0; { - int prevk = -1,prevl = -1; - for (auto [k,l,v] : subklv) { + int prevk = -1, prevl = -1; + for (auto [k, l, v] : subklv) { if (prevk != k || prevl != l) { - ++nunique; prevk = k; prevl = l; + ++nunique; + prevk = k; + prevl = l; } } } - - subk.clear(); subk.reserve(nunique); - subl.clear(); subl.reserve(nunique); - val.clear(); val.reserve(nunique); - int prevk = -1,prevl = -1; - for (auto [k,l,v] : subklv) { + subk.clear(); + subk.reserve(nunique); + subl.clear(); + subl.reserve(nunique); + val.clear(); + val.reserve(nunique); + + int prevk = -1, prevl = -1; + for (auto [k, l, v] : subklv) { if (prevk == k && prevl == l) { val.back() += v; - } - else { - subk.push_back(k); prevk = k; - subl.push_back(l); prevl = l; + } else { + subk.push_back(k); + prevk = k; + subl.push_back(l); + prevl = l; val.push_back(v); } } } } -absl::Status MosekSolver::AddConstraint(int64_t id, const QuadraticConstraintProto& cons) { +absl::Status MosekSolver::AddConstraint(int64_t id, + const QuadraticConstraintProto& cons) { int coni = msk.NumCon(); double clb = cons.lower_bound(); double cub = cons.upper_bound(); @@ -294,12 +292,15 @@ absl::Status MosekSolver::AddConstraint(int64_t id, const QuadraticConstraintPro } size_t nnz = cons.linear_terms().ids_size(); - std::vector subj; subj.reserve(nnz); - std::vector val; val.reserve(nnz); + std::vector subj; + subj.reserve(nnz); + std::vector val; + val.reserve(nnz); - for (const auto id : cons.linear_terms().ids()) subj.push_back(variable_map[id]); + for (const auto id : cons.linear_terms().ids()) + subj.push_back(variable_map[id]); for (const auto c : cons.linear_terms().values()) val.push_back(c); - RETURN_IF_ERROR(msk.PutARow(coni,subj, val)); + RETURN_IF_ERROR(msk.PutARow(coni, subj, val)); // quadratic terms @@ -308,9 +309,9 @@ absl::Status MosekSolver::AddConstraint(int64_t id, const QuadraticConstraintPro std::vector subl; std::vector val; - SparseDoubleMatrixToTril(cons.quadratic_terms(),subk,subl,val); + SparseDoubleMatrixToTril(cons.quadratic_terms(), subk, subl, val); - RETURN_IF_ERROR(msk.PutQCon(coni,subk,subl,val)); + RETURN_IF_ERROR(msk.PutQCon(coni, subk, subl, val)); } quadconstr_map[id] = coni; @@ -466,10 +467,8 @@ absl::Status MosekSolver::AddConicConstraints( } sizes.push_back(expr0.ids_size()); - for (const auto& id : expr0.ids()) - subj.push_back(variable_map[id]); - for (auto c : expr0.coefficients()) - cof.push_back(c); + for (const auto& id : expr0.ids()) subj.push_back(variable_map[id]); + for (auto c : expr0.coefficients()) cof.push_back(c); for (const auto& expri : con.arguments_to_norm()) { sizes.push_back(expri.ids_size()); @@ -524,15 +523,17 @@ absl::StatusOr MosekSolver::Update(const ModelUpdateProto& model_update) { } } - for (auto id : model_update.quadratic_constraint_updates().deleted_constraint_ids()) { + for (auto id : + model_update.quadratic_constraint_updates().deleted_constraint_ids()) { if (quadconstr_map.contains(id)) { int i = quadconstr_map[id]; quadconstr_map.erase(id); RETURN_IF_ERROR(msk.ClearConstraint(i)); } } - - for (auto &[id,con] : model_update.quadratic_constraint_updates().new_constraints()) { + + for (auto& [id, con] : + model_update.quadratic_constraint_updates().new_constraints()) { RETURN_IF_ERROR(AddConstraint(id, con)); } @@ -617,19 +618,17 @@ absl::Status MosekSolver::UpdateObjective( for (auto id : objupds.linear_coefficients().ids()) subj.push_back(variable_map[id]); - if (objupds.quadratic_coefficients().column_ids_size() > 0) { // note: this specifies the full q matrix, and we assume that it is // symmetric and only specifies the lower triangular part. - auto & qobj = objupds.quadratic_coefficients(); + auto& qobj = objupds.quadratic_coefficients(); size_t nqnz = qobj.row_ids_size(); - std::vector> subklv; subklv.reserve(nqnz); - for (auto [kit,lit,cit] = std::make_tuple( - qobj.row_ids().begin(), - qobj.column_ids().begin(), - qobj.coefficients().begin()); - kit != qobj.row_ids().end() && - lit != qobj.column_ids().end() && + std::vector> subklv; + subklv.reserve(nqnz); + for (auto [kit, lit, cit] = + std::make_tuple(qobj.row_ids().begin(), qobj.column_ids().begin(), + qobj.coefficients().begin()); + kit != qobj.row_ids().end() && lit != qobj.column_ids().end() && cit != qobj.coefficients().end(); ++kit, ++lit, ++cit) { int k = variable_map[*kit]; @@ -637,38 +636,35 @@ absl::Status MosekSolver::UpdateObjective( int v = variable_map[*cit]; if (k < l) { - subklv.push_back({l,k,v}); - } - else { - subklv.push_back({k,l,v}); + subklv.push_back({l, k, v}); + } else { + subklv.push_back({k, l, v}); } } - std::sort(subklv.begin(),subklv.end()); - - std::vector subk; subk.reserve(nqnz); - std::vector subl; subl.reserve(nqnz); - std::vector val; val.reserve(nqnz); + std::sort(subklv.begin(), subklv.end()); + + std::vector subk; + subk.reserve(nqnz); + std::vector subl; + subl.reserve(nqnz); + std::vector val; + val.reserve(nqnz); - int prevk = -1,prevl = -1; - for (auto [k,l,v] : subklv) { + int prevk = -1, prevl = -1; + for (auto [k, l, v] : subklv) { if (prevk == k && prevl == l) { val.back() += v; - } - else { - subk.push_back(k); prevk = k; - subk.push_back(l); prevl = l; + } else { + subk.push_back(k); + prevk = k; + subk.push_back(l); + prevl = l; val.push_back(v); } } - } - - - - - RETURN_IF_ERROR(msk.UpdateObjectiveSense(objupds.direction_update())); RETURN_IF_ERROR(msk.UpdateObjective(objupds.offset_update(), subj, cof)); @@ -716,15 +712,16 @@ absl::StatusOr> MosekSolver::New( RETURN_IF_ERROR(mskslv->ReplaceObjective(model.objective())); RETURN_IF_ERROR(mskslv->AddConstraints(model.linear_constraints(), model.linear_constraint_matrix())); - RETURN_IF_ERROR(mskslv->AddConicConstraints(model.second_order_cone_constraints())); - for (auto &[k,v] : model.quadratic_constraints()) - RETURN_IF_ERROR(mskslv->AddConstraint(k,v)); + RETURN_IF_ERROR( + mskslv->AddConicConstraints(model.second_order_cone_constraints())); + for (auto& [k, v] : model.quadratic_constraints()) + RETURN_IF_ERROR(mskslv->AddConstraint(k, v)); RETURN_IF_ERROR( mskslv->AddIndicatorConstraints(model.indicator_constraints())); std::unique_ptr res(std::move(mskslv)); - //std::cout << " ---------- MosekSolver::New" << std::endl; + // std::cout << " ---------- MosekSolver::New" << std::endl; return res; } @@ -1000,7 +997,7 @@ absl::StatusOr MosekSolver::Solve( // - EmphasisProto cuts // - EmphasisProto heuristics // - EmphasisProto scaling - + auto solve_start = absl::Now(); double dpar_optimizer_max_time = msk.GetParam(MSK_DPAR_OPTIMIZER_MAX_TIME); @@ -1154,7 +1151,7 @@ absl::StatusOr MosekSolver::Solve( { BufferedMessageCallback bmsg_cb(message_cb); // TODO: Use model_parameters - //msk.WriteData("test.ptf"); + // msk.WriteData("test.ptf"); auto r = msk.Optimize( [&](const std::string& msg) { bmsg_cb.OnMessage(msg); }, [&](MSKcallbackcodee code, absl::Span dinf, @@ -1267,13 +1264,13 @@ absl::StatusOr MosekSolver::Solve( if (solsta == mosek::SolSta::OPTIMAL || solsta == mosek::SolSta::INTEGER_OPTIMAL) { - trmp = OptimalTerminationProto(msk.GetPrimalObj(whichsol), - msk.GetDualObj(whichsol), - ""); - // Hack: - trmp.mutable_objective_bounds()->set_primal_bound(trmp.objective_bounds().primal_bound()); - trmp.mutable_objective_bounds()->set_dual_bound(trmp.objective_bounds().dual_bound()); + msk.GetDualObj(whichsol), ""); + // Hack: + trmp.mutable_objective_bounds()->set_primal_bound( + trmp.objective_bounds().primal_bound()); + trmp.mutable_objective_bounds()->set_dual_bound( + trmp.objective_bounds().dual_bound()); } else if (solsta == mosek::SolSta::PRIM_INFEAS_CER) trmp = InfeasibleTerminationProto( @@ -1341,16 +1338,14 @@ absl::StatusOr MosekSolver::Solve( case mosek::SolSta::INTEGER_OPTIMAL: case mosek::SolSta::PRIM_FEAS: case mosek::SolSta::DUAL_FEAS: - case mosek::SolSta::PRIM_AND_DUAL_FEAS: - { - auto r = Solution(whichsol, ordered_xc_ids, ordered_xx_ids, - skip_xx_zeros, ordered_y_ids, skip_y_zeros, - ordered_yx_ids, skip_yx_zeros); - if (r.ok()) { - *result.add_solutions() = std::move(*r); - } + case mosek::SolSta::PRIM_AND_DUAL_FEAS: { + auto r = Solution(whichsol, ordered_xc_ids, ordered_xx_ids, + skip_xx_zeros, ordered_y_ids, skip_y_zeros, + ordered_yx_ids, skip_yx_zeros); + if (r.ok()) { + *result.add_solutions() = std::move(*r); } - break; + } break; case mosek::SolSta::DUAL_INFEAS_CER: { auto r = PrimalRay(whichsol, ordered_xx_ids, skip_xx_zeros); if (r.ok()) { @@ -1371,12 +1366,14 @@ absl::StatusOr MosekSolver::Solve( } } - SolveStatsProto * stats = result.mutable_solve_stats(); - stats->set_simplex_iterations(msk.GetIntInfoItem(MSK_IINF_SIM_PRIMAL_ITER)+msk.GetIntInfoItem(MSK_IINF_SIM_DUAL_ITER)); + SolveStatsProto* stats = result.mutable_solve_stats(); + stats->set_simplex_iterations(msk.GetIntInfoItem(MSK_IINF_SIM_PRIMAL_ITER) + + msk.GetIntInfoItem(MSK_IINF_SIM_DUAL_ITER)); stats->set_barrier_iterations(msk.GetIntInfoItem(MSK_IINF_INTPNT_ITER)); stats->set_node_count(msk.GetIntInfoItem(MSK_IINF_MIO_NUM_SOLVED_NODES)); - auto solve_time_proto = util_time::EncodeGoogleApiProto(absl::Now() - solve_start); + auto solve_time_proto = + util_time::EncodeGoogleApiProto(absl::Now() - solve_start); if (solve_time_proto.ok()) { *stats->mutable_solve_time() = *solve_time_proto; } @@ -1394,4 +1391,4 @@ MosekSolver::ComputeInfeasibleSubsystem(const SolveParametersProto&, MATH_OPT_REGISTER_SOLVER(SOLVER_TYPE_MOSEK, MosekSolver::New); -} // namespace operations_research::math_opt +} // namespace operations_research::math_opt diff --git a/ortools/math_opt/solvers/mosek_solver.h b/ortools/math_opt/solvers/mosek_solver.h index ff2b2afa1f5..14e012c257c 100644 --- a/ortools/math_opt/solvers/mosek_solver.h +++ b/ortools/math_opt/solvers/mosek_solver.h @@ -103,10 +103,9 @@ class MosekSolver : public SolverInterface { bool skip_y_zeros, const std::vector& ordered_yx_ids, bool skip_yx_zeros); - void SparseDoubleMatrixToTril(const SparseDoubleMatrixProto & qdata, - std::vector & subi, - std::vector & subj, - std::vector & cof); + void SparseDoubleMatrixToTril(const SparseDoubleMatrixProto& qdata, + std::vector& subi, std::vector& subj, + std::vector& cof); MosekSolver(Mosek&& msk); MosekSolver(MosekSolver&) = delete; From 8e3f9db1115a3ad2c5f0794e5d17e8964d836627 Mon Sep 17 00:00:00 2001 From: Ulfw Date: Wed, 30 Oct 2024 17:11:26 +0100 Subject: [PATCH 09/10] Update copyright notices --- ortools/linear_solver/mosek_interface.cc | 2 +- ortools/math_opt/solvers/mosek/mosekwrp.cc | 15 +++++++++++++++ ortools/math_opt/solvers/mosek/mosekwrp.h | 15 +++++++++++++++ ortools/math_opt/solvers/mosek_solver.cc | 2 +- ortools/math_opt/solvers/mosek_solver.h | 2 +- 5 files changed, 33 insertions(+), 3 deletions(-) diff --git a/ortools/linear_solver/mosek_interface.cc b/ortools/linear_solver/mosek_interface.cc index 5c69c0496ff..b3726f6fc1d 100644 --- a/ortools/linear_solver/mosek_interface.cc +++ b/ortools/linear_solver/mosek_interface.cc @@ -1,4 +1,4 @@ -// Copyright 2010-2024 Google LLg +// Copyright 2024 Mosek // 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 diff --git a/ortools/math_opt/solvers/mosek/mosekwrp.cc b/ortools/math_opt/solvers/mosek/mosekwrp.cc index 85e51a2f48e..4cce003fa8b 100644 --- a/ortools/math_opt/solvers/mosek/mosekwrp.cc +++ b/ortools/math_opt/solvers/mosek/mosekwrp.cc @@ -1,3 +1,18 @@ +// Copyright 2024 Mosek +// 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. + +// Mosek backend to MPSolver. +// #include "mosekwrp.h" #include diff --git a/ortools/math_opt/solvers/mosek/mosekwrp.h b/ortools/math_opt/solvers/mosek/mosekwrp.h index 70fb5f7d644..5a2959bb60a 100644 --- a/ortools/math_opt/solvers/mosek/mosekwrp.h +++ b/ortools/math_opt/solvers/mosek/mosekwrp.h @@ -1,3 +1,18 @@ +// Copyright 2024 Mosek +// 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. + +// Mosek backend to MPSolver. +// #ifndef OR_TOOLS_MATH_OPT_SOLVERS_MOSEK_MOSEKWRP_H_ #define OR_TOOLS_MATH_OPT_SOLVERS_MOSEK_MOSEKWRP_H_ diff --git a/ortools/math_opt/solvers/mosek_solver.cc b/ortools/math_opt/solvers/mosek_solver.cc index 589d1df56b2..28d92939f06 100644 --- a/ortools/math_opt/solvers/mosek_solver.cc +++ b/ortools/math_opt/solvers/mosek_solver.cc @@ -1,4 +1,4 @@ -// Copyright 2010-2024 Mosek +// Copyright 2024 Mosek // 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 diff --git a/ortools/math_opt/solvers/mosek_solver.h b/ortools/math_opt/solvers/mosek_solver.h index 14e012c257c..2e864825f20 100644 --- a/ortools/math_opt/solvers/mosek_solver.h +++ b/ortools/math_opt/solvers/mosek_solver.h @@ -1,4 +1,4 @@ -// Copyright 2010-2024 Google LLC +// Copyright 2024 Mosek // 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 From 16fc39841fc863636a3becdeb73a155a7c7ed636 Mon Sep 17 00:00:00 2001 From: ulfw Date: Fri, 1 Nov 2024 11:27:44 +0100 Subject: [PATCH 10/10] fix mosek/cmake config for windows --- cmake/FindMOSEK.cmake | 48 ++++++++++++++++++++++++++++++++++--------- 1 file changed, 38 insertions(+), 10 deletions(-) diff --git a/cmake/FindMOSEK.cmake b/cmake/FindMOSEK.cmake index 515c292266a..3e2601b2fef 100644 --- a/cmake/FindMOSEK.cmake +++ b/cmake/FindMOSEK.cmake @@ -100,6 +100,26 @@ function(FindMosekPlatformInPath RESULT PATH) endfunction() +function(MosekVersionFromPath RESVMAJOR RESVMINOR PATH) + execute_process(COMMAND "${PATH}/bin/mosek" "-v" RESULT_VARIABLE CMDRES OUTPUT_VARIABLE VERTEXT) + if (${CMDRES} EQUAL 0) + if (${VERTEXT} MATCHES "^MOSEK version [0-9]+[.][0-9]+[.][0-9]+") + string(REGEX REPLACE "^MOSEK version ([0-9]+)[.]([0-9]+).*" "\\1;\\2" MSKVER ${VERTEXT}) + + list(GET MSKVER 0 VMAJOR) + list(GET MSKVER 1 VMINOR) + set(${RESVMAJOR} "${VMAJOR}" PARENT_SCOPE) + set(${RESVMINOR} "${VMINOR}" PARENT_SCOPE) + + #message(STATUS "mskver = ${VMAJOR}.${VMINOR} -> ${${RESVMAJOR}}.${${RESVMINOR}}") + #message(STATUS " ${RESVMAJOR} = ${${RESVMAJOR}}") + #message(STATUS " ${RESVMINOR} = ${${RESVMINOR}}") + return() + endif() + endif() +endfunction() + + # Where to look for MOSEK: # Standard procedure in Linux/OSX is to install MOSEK in the home directory, i.e. # $HOME/mosek/X.Y/... @@ -129,24 +149,32 @@ if(NOT MOSEK_PLATFORM_DIR) endif() if(NOT MOSEK_PLATFORM_DIR) # Look in users home dir - if(ENV{HOMEDRIVE} AND ENV{HOMEPATH}) - FindMosekPlatformInPath(MOSEK_PLATFORM_DIR"$ENV{HOMEDRIVE}$ENV{HOMEPATH}") + if(EXISTS "$ENV{HOMEDRIVE}$ENV{HOMEPATH}") + FindMosekPlatformInPath(MOSEK_PLATFORM_DIR "$ENV{HOMEDRIVE}$ENV{HOMEPATH}") endif() endif() if(NOT MOSEK_PLATFORM_DIR) - message(FATAL_ERROR " MOSEK_PLATFORM_DIR could not be detected") + message(FATAL_ERROR " MOSEK_PLATFORM_DIR could not be detected") else() - message(STATUS " MOSEK_PLATFORM_DIR detected: ${MOSEK_PLATFORM_DIR}") - set(MOSEK_PFDIR_FOUND TRUE) + MosekVersionFromPath(MSKVMAJOR MSKVMINOR "${MOSEK_PLATFORM_DIR}") + if(MSKVMAJOR) + message(STATUS " MOSEK ${MSKVMAJOR}.${MSKVMINOR}: ${MOSEK_PLATFORM_DIR}") + set(MOSEK_PFDIR_FOUND TRUE) + else() + message(STATUS " Found MOSEK_PLATFORM_DIR, but failed to determine version") + endif() endif() - -if(MOSEK_PFDIR_FOUND AND NOT TARGET Mosek::Mosek) +if(MOSEK_PFDIR_FOUND AND MSKVMAJOR AND MSKVMINOR AND NOT TARGET mosek::mosek) add_library(mosek::mosek UNKNOWN IMPORTED) find_path(MOSEKINC mosek.h HINTS ${MOSEK_PLATFORM_DIR}/h) - find_library(LIBMOSEK mosek64 HINTS ${MOSEK_PLATFORM_DIR}/bin) + if(WIN32) + find_library(LIBMOSEK mosek64_ HINTS ${MOSEK_PLATFORM_DIR}/bin) + else() + find_library(LIBMOSEK mosek64 HINTS ${MOSEK_PLATFORM_DIR}/bin) + endif() if(LIBMOSEK) set_target_properties(mosek::mosek PROPERTIES IMPORTED_LOCATION ${LIBMOSEK}) @@ -155,7 +183,7 @@ if(MOSEK_PFDIR_FOUND AND NOT TARGET Mosek::Mosek) if (MOSEKINC) target_include_directories(mosek::mosek INTERFACE "${MOSEKINC}") endif() -elseif(NOT TARGET Mosek::Mosek) +elseif(NOT TARGET mosek::mosek) add_library(mosek::mosek UNKNOWN IMPORTED) find_path(MOSEKINC mosek.h) @@ -166,6 +194,6 @@ elseif(NOT TARGET Mosek::Mosek) endif() if (MOSEKINC) - target_include_directories(mosek::mosek INTERFACE "${MOSEKINC}") + target_include_directories(mosek::mosek INTERFACE "${MOSEKINC}") endif() endif()