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..3e2601b2fef --- /dev/null +++ b/cmake/FindMOSEK.cmake @@ -0,0 +1,199 @@ +# 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() + + +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/... +# 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(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") +else() + 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 MSKVMAJOR AND MSKVMINOR AND NOT TARGET mosek::mosek) + add_library(mosek::mosek UNKNOWN IMPORTED) + + find_path(MOSEKINC mosek.h HINTS ${MOSEK_PLATFORM_DIR}/h) + 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}) + 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..43b3833d2e0 100644 --- a/cmake/check_deps.cmake +++ b/cmake/check_deps.cmake @@ -116,6 +116,14 @@ if(USE_CPLEX AND NOT TARGET CPLEX::CPLEX) message(FATAL_ERROR "Target CPLEX::CPLEX not available.") endif() +# Check optional Dependencies +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 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/portfolio_1_basic.cc b/examples/cpp/portfolio_1_basic.cc new file mode 100644 index 00000000000..3b3df9919c3 --- /dev/null +++ b/examples/cpp/portfolio_1_basic.cc @@ -0,0 +1,101 @@ +#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::SolverType st) { + + 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)); + + 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) { + 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, st, 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,math_opt::SolverType::kMosek); + + if (! res.ok()) { + std::cerr << "Failed to solve problem: " << res.status() << 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; + } +} 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 2e13e952b6e..dd3c9df6488 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..b3726f6fc1d --- /dev/null +++ b/ortools/linear_solver/mosek_interface.cc @@ -0,0 +1,1402 @@ +// 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 + +#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 bdfa1470de1..4deba851e6d 100644 --- a/ortools/math_opt/solvers/BUILD.bazel +++ b/ortools/math_opt/solvers/BUILD.bazel @@ -103,6 +103,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..4cce003fa8b --- /dev/null +++ b/ortools/math_opt/solvers/mosek/mosekwrp.cc @@ -0,0 +1,686 @@ +// 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 + +#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"); + + 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] + ? 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( + 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"); + 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 (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] + ? 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, ub); + 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::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())) { + 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, + 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) { + auto [msg,rname,r] = std::move(LastError()); + return absl::InternalError((std::stringstream() << rname << "(" << r << "): " << msg).str()); + } + 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::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())) { + auto [msg,rname,r] = std::move(LastError()); + return absl::InternalError((std::stringstream() << rname << "(" << r << "): " << msg).str()); + } + + 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())) { + 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) { + 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())) { + auto [msg,rname,r] = std::move(LastError()); + return absl::InternalError((std::stringstream() << rname << "(" << r << "): " << msg).str()); + } + 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()) { + 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: + 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); + MSK_putqconk(task.get(), i, 0, &subj,&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) { + std::cerr << "MOSEK ERROR: " << r << std::endl; + return absl::InternalError("Optimization failedcode"); + } + + 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) { + std::cerr << "MOSEK ERROR: " << r << std::endl; + 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()); +} + +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); +} +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; +} + +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); + 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..5a2959bb60a --- /dev/null +++ b/ortools/math_opt/solvers/mosek/mosekwrp.h @@ -0,0 +1,217 @@ +// 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_ + +#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(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); + + 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; + + void GetC(std::vector & c, double & cfix); + + 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; + } + + + int GetIntInfoItem(MSKiinfiteme item); + int64_t GetLongInfoItem(MSKliinfiteme item); + double GetDoubleInfoItem(MSKdinfiteme item); + + 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); + + absl::StatusOr append_cons(int n, const double* lb, + const double* ub); +}; + +} // 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..28d92939f06 --- /dev/null +++ b/ortools/math_opt/solvers/mosek_solver.cc @@ -0,0 +1,1394 @@ +// 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. + +#include "mosek_solver.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#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/math_opt_proto_utils.h" +#include "ortools/math_opt/core/solver_interface.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/math_opt/sparse_containers.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, + .quadratic_objectives = SupportType::kSupported, + .quadratic_constraints = 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)); + + // 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)); + } + + 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]; + double v = *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; + subl.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(); + 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) { + 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()) { + totalnnz += lexp.ids_size(); + } + + subj.reserve(totalnnz); + cof.reserve(totalnnz); + + b.push_back(expr0.offset()); + for (const auto& expri : con.arguments_to_norm()) { + b.push_back(expri.offset()); + } + + 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& 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); + } + } + + 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()) { + 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()) { + 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()) { + 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()) { + 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())); + 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]); + + 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)); + + 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.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->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; + + 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.GetDualObj(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 + + 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); + 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 = MSK_RES_OK; + { + BufferedMessageCallback bmsg_cb(message_cb); + // TODO: Use model_parameters + // msk.WriteData("test.ptf"); + 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 { + 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) { + 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()); + + } 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; + } + } + + 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; +} + +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..2e864825f20 --- /dev/null +++ b/ortools/math_opt/solvers/mosek_solver.h @@ -0,0 +1,116 @@ +// 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. + +#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 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); + 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); + + void SparseDoubleMatrixToTril(const SparseDoubleMatrixProto& qdata, + std::vector& subi, std::vector& subj, + std::vector& cof); + + MosekSolver(Mosek&& msk); + MosekSolver(MosekSolver&) = delete; +}; + +} // 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.