From f50fc893a8c3ac616b54ff2e374ede1c183a2ed9 Mon Sep 17 00:00:00 2001 From: Manuel Ruiz Date: Tue, 6 Aug 2019 18:21:29 +0200 Subject: [PATCH] xpress integration in or-tools --- CMakeLists.txt | 79 + cmake/dependencies/CMakeLists.txt | 1 + ...trawberry_fields_with_column_generation.cc | 15 + ortools/algorithms/knapsack_solver.cc | 12 + ortools/algorithms/knapsack_solver.h | 7 +- ortools/linear_solver/CMakeLists.txt | 20 + ortools/linear_solver/linear_solver.cc | 19 + ortools/linear_solver/linear_solver.h | 10 + ortools/linear_solver/xpress_interface.cc | 1485 +++++++++++++++++ 9 files changed, 1647 insertions(+), 1 deletion(-) create mode 100644 ortools/linear_solver/xpress_interface.cc diff --git a/CMakeLists.txt b/CMakeLists.txt index 29357c1dc6..66efdaeca2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,6 +7,14 @@ if(POLICY CMP0068) cmake_policy(SET CMP0068 NEW) endif() +if(POLICY CMP0078) + cmake_policy(SET CMP0078 OLD) +endif() + +if(POLICY CMP0086) + cmake_policy(SET CMP0086 NEW) +endif() + include(utils) set_version(VERSION) @@ -32,10 +40,14 @@ option(BUILD_CXX "Build C++ library" ON) option(BUILD_PYTHON "Build Python Library" OFF) option(BUILD_JAVA "Build Java Library" OFF) option(BUILD_DOTNET "Build .NET Library" OFF) +option(USE_XPRESS "Build and use XPRESS interface" OFF) +option(USE_CPLEX "Build and use CPLEX interface" OFF) message(STATUS "Build C++ library: ${BUILD_CXX}") message(STATUS "Build Python: ${BUILD_PYTHON}") message(STATUS "Build Java: ${BUILD_JAVA}") message(STATUS "Build .Net: ${BUILD_DOTNET}") +message(STATUS "USE_XPRESS: ${USE_XPRESS}") +message(STATUS "USE_CPLEX: ${USE_CPLEX}") #option(BUILD_DOC "Build doxygen" OFF) #message(STATUS "Build doxygen: ${BUILD_DOC}") @@ -79,6 +91,49 @@ message(STATUS "Build Clp: ${BUILD_Clp}") message(STATUS "Build Cgl: ${BUILD_Cgl}") message(STATUS "Build Cbc: ${BUILD_Cbc}") +if (USE_XPRESS) + + if (APPLE) + message(FATAL_ERROR "XPRESS not yet supported on MACOS") + endif() + + message(STATUS "XPRESS: configuring makefiles") + if (NOT XPRESSDIR) + set(XPRESSDIR $ENV{XPRESSDIR}) + endif(NOT XPRESSDIR) + + message(STATUS "XPRESSDIR: ${XPRESSDIR}") + + if (NOT XPRESSDIR) + message(FATAL_ERROR "XPRESSDIR: not found") + endif() + + add_compile_definitions(USE_XPRESS) +endif(USE_XPRESS) + + +if (USE_CPLEX) + + if (APPLE) + message(FATAL_ERROR "CPLEX not yet supported on MACOS") + elseif(UNIX) + message(FATAL_ERROR "CPLEX not yet supported on LINUX") + endif() + + message(STATUS "CPLEX: configuring makefiles") + if (NOT CPLEXDIR) + set(CPLEXDIR $ENV{CPLEXDIR}) + endif(NOT CPLEXDIR) + + message(STATUS "CPLEXDIR: ${CPLEXDIR}") + + if (NOT CPLEXDIR) + message(FATAL_ERROR "CPLEXDIR: not found") + endif() + + add_compile_definitions(USE_CPLEX) +endif(USE_CPLEX) + # Build Needed dependencies add_subdirectory(cmake/dependencies dependencies) list(APPEND CMAKE_PREFIX_PATH ${CMAKE_BINARY_DIR}/dependencies/install) @@ -90,3 +145,27 @@ include(cpp) include(python) include(java) include(dotnet) + +if (USE_XPRESS) + if (APPLE) + message(FATAL_ERROR "XPRESS not yet supported on MACOS") + #target_link_libraries(ortools PUBLIC ${XPRESSDIR}/lib/libxprs.dylib) + elseif(UNIX) + target_link_libraries(ortools PUBLIC ${XPRESSDIR}/lib/libxprs.so) + elseif(MSVC) + target_link_libraries(ortools PUBLIC ${XPRESSDIR}/lib/xprs.lib) + endif() +endif(USE_XPRESS) + + +if (USE_CPLEX) + if (APPLE) + message(FATAL_ERROR "CPLEX not yet supported on MACOS") + #target_link_libraries(ortools PUBLIC ${CPLEXDIR}/lib/x64_windows_vs2015/stat_mda/cplex1270.dylib) + elseif(UNIX) + message(FATAL_ERROR "CPLEX not yet supported on LINUX") + #target_link_libraries(ortools PUBLIC ${CPLEXDIR}/lib/x64_windows_vs2015/stat_mda/cplex1270.so) + elseif(MSVC) + target_link_libraries(ortools PUBLIC ${CPLEXDIR}/lib/x64_windows_vs2015/stat_mda/cplex1270.lib) + endif() +endif(USE_CPLEX) diff --git a/cmake/dependencies/CMakeLists.txt b/cmake/dependencies/CMakeLists.txt index d6ad4c1857..d6c29cebf6 100644 --- a/cmake/dependencies/CMakeLists.txt +++ b/cmake/dependencies/CMakeLists.txt @@ -104,6 +104,7 @@ if(BUILD_Protobuf) TAG "v3.8.0" CMAKE_ARGS + -Dprotobuf_MSVC_STATIC_RUNTIME:BOOL=OFF "SOURCE_SUBDIR cmake" ) endif() diff --git a/examples/cpp/strawberry_fields_with_column_generation.cc b/examples/cpp/strawberry_fields_with_column_generation.cc index 77aba824d4..024f58cc5a 100644 --- a/examples/cpp/strawberry_fields_with_column_generation.cc +++ b/examples/cpp/strawberry_fields_with_column_generation.cc @@ -618,11 +618,26 @@ int main(int argc, char** argv) { found = true; } #endif // USE_GLOP +#if defined(USE_XPRESS) + if (FLAGS_colgen_solver == "xpress") { + solver_type = operations_research::MPSolver::XPRESS_LINEAR_PROGRAMMING; + //solver_type = operations_research::MPSolver::CPLEX_LINEAR_PROGRAMMING; + found = true; + } +#endif +#if defined(USE_CPLEX) + if (FLAGS_colgen_solver == "cplex") { + solver_type = operations_research::MPSolver::CPLEX_LINEAR_PROGRAMMING; + found = true; + } +#endif if (!found) { LOG(ERROR) << "Unknown solver " << FLAGS_colgen_solver; return 1; } + LOG(INFO) << "Chosen solver: " << FLAGS_colgen_solver << std::endl; + if (FLAGS_colgen_instance == -1) { for (int i = 0; i < operations_research::kInstanceCount; ++i) { const operations_research::Instance& instance = diff --git a/ortools/algorithms/knapsack_solver.cc b/ortools/algorithms/knapsack_solver.cc index d8a2a08ae0..6c1b6e1577 100644 --- a/ortools/algorithms/knapsack_solver.cc +++ b/ortools/algorithms/knapsack_solver.cc @@ -1153,6 +1153,18 @@ KnapsackSolver::KnapsackSolver(SolverType solver_type, MPSolver::SCIP_MIXED_INTEGER_PROGRAMMING, solver_name); break; #endif // USE_SCIP +#if defined(USE_XPRESS) + case KNAPSACK_MULTIDIMENSION_XPRESS_MIP_SOLVER: + solver_ = absl::make_unique( + MPSolver::XPRESS_MIXED_INTEGER_PROGRAMMING, solver_name); + break; +#endif +#if defined(USE_CPLEX) + case KNAPSACK_MULTIDIMENSION_CPLEX_MIP_SOLVER: + solver_ = absl::make_unique( + MPSolver::CPLEX_MIXED_INTEGER_PROGRAMMING, solver_name); + break; +#endif default: LOG(FATAL) << "Unknown knapsack solver type."; } diff --git a/ortools/algorithms/knapsack_solver.h b/ortools/algorithms/knapsack_solver.h index be53e4ac30..2639d170d9 100644 --- a/ortools/algorithms/knapsack_solver.h +++ b/ortools/algorithms/knapsack_solver.h @@ -157,7 +157,6 @@ class KnapsackSolver { */ KNAPSACK_MULTIDIMENSION_CBC_MIP_SOLVER = 3, #endif // USE_CBC - /** Generic Solver. * * This solver can deal with both large number of items and several @@ -173,6 +172,12 @@ class KnapsackSolver { */ KNAPSACK_MULTIDIMENSION_SCIP_MIP_SOLVER = 6, #endif // USE_SCIP +#if defined(USE_XPRESS) + KNAPSACK_MULTIDIMENSION_XPRESS_MIP_SOLVER = 7, +#endif +#if defined(USE_CPLEX) + KNAPSACK_MULTIDIMENSION_CPLEX_MIP_SOLVER = 8, +#endif }; explicit KnapsackSolver(const std::string& solver_name); diff --git a/ortools/linear_solver/CMakeLists.txt b/ortools/linear_solver/CMakeLists.txt index bfdb78b8eb..35f8239fcc 100644 --- a/ortools/linear_solver/CMakeLists.txt +++ b/ortools/linear_solver/CMakeLists.txt @@ -58,3 +58,23 @@ add_dependencies(${NAME} Coin::Cbc ${PROJECT_NAME}::proto) add_library(${PROJECT_NAME}::linear_solver ALIAS ${NAME}) + +if (USE_XPRESS) + target_include_directories(${NAME} PUBLIC ${XPRESSDIR}/include) + if(UNIX) + target_link_libraries(${NAME} PUBLIC ${XPRESSDIR}/lib/libxprs.so) + elseif(MSVC) + target_link_libraries(${NAME} PUBLIC ${XPRESSDIR}/lib/xprs.lib) + endif() +endif(USE_XPRESS) + +if (USE_CPLEX) + target_include_directories(${NAME} PUBLIC + ${CPLEXDIR}/include + ) + if(UNIX) + target_link_libraries(${NAME} PUBLIC ${CPLEXDIR}/lib/x64_windows_vs2015/stat_mda/cplex1270.a) + elseif(MSVC) + target_link_libraries(${NAME} PUBLIC ${CPLEXDIR}/lib/x64_windows_vs2015/stat_mda/cplex1270.lib) + endif() +endif(USE_CPLEX) \ No newline at end of file diff --git a/ortools/linear_solver/linear_solver.cc b/ortools/linear_solver/linear_solver.cc index 9aa5727224..bb2b407481 100644 --- a/ortools/linear_solver/linear_solver.cc +++ b/ortools/linear_solver/linear_solver.cc @@ -349,6 +349,9 @@ extern MPSolverInterface* BuildCplexInterface(bool mip, MPSolver* const solver); extern MPSolverInterface* BuildGLOPInterface(MPSolver* const solver); #endif +#if defined(USE_XPRESS) +extern MPSolverInterface* BuildXpressInterface(bool mip, MPSolver* const solver); +#endif namespace { MPSolverInterface* BuildSolverInterface(MPSolver* const solver) { @@ -391,6 +394,12 @@ MPSolverInterface* BuildSolverInterface(MPSolver* const solver) { return BuildCplexInterface(false, solver); case MPSolver::CPLEX_MIXED_INTEGER_PROGRAMMING: return BuildCplexInterface(true, solver); +#endif +#if defined(USE_XPRESS) + case MPSolver::XPRESS_MIXED_INTEGER_PROGRAMMING: + return BuildXpressInterface(true, solver); + case MPSolver::XPRESS_LINEAR_PROGRAMMING: + return BuildXpressInterface(false, solver); #endif default: // TODO(user): Revert to the best *available* interface. @@ -455,6 +464,14 @@ bool MPSolver::SupportsProblemType(OptimizationProblemType problem_type) { #endif #ifdef USE_CBC if (problem_type == CBC_MIXED_INTEGER_PROGRAMMING) return true; +#endif +#ifdef USE_XPRESS + if (problem_type == XPRESS_MIXED_INTEGER_PROGRAMMING) return true; + if (problem_type == XPRESS_LINEAR_PROGRAMMING) return true; +#endif +#ifdef USE_CPLEX + if (problem_type == CPLEX_LINEAR_PROGRAMMING) return true; + if (problem_type == CPLEX_MIXED_INTEGER_PROGRAMMING) return true; #endif return false; } @@ -1458,7 +1475,9 @@ void MPSolver::GenerateConstraintNameIndex() const { } } +#if defined(USE_GUROBI) bool MPSolver::NextSolution() { return interface_->NextSolution(); } +#endif // ---------- MPSolverInterface ---------- diff --git a/ortools/linear_solver/linear_solver.h b/ortools/linear_solver/linear_solver.h index 13fbd253fd..9413c79a7b 100644 --- a/ortools/linear_solver/linear_solver.h +++ b/ortools/linear_solver/linear_solver.h @@ -224,6 +224,10 @@ class MPSolver { #if defined(USE_BOP) /// Linear Boolean Programming Solver. BOP_INTEGER_PROGRAMMING = 12, +#endif +#if defined(USE_XPRESS) + XPRESS_LINEAR_PROGRAMMING = 101, + XPRESS_MIXED_INTEGER_PROGRAMMING = 102, #endif }; @@ -726,8 +730,10 @@ class MPSolver { * As of 2018-08-09, only Gurobi supports NextSolution(), see * linear_solver_underlying_gurobi_test for an example of how to configure * Gurobi for this purpose. The other solvers return false unconditionally. +#if defined(USE_GUROBI) */ ABSL_MUST_USE_RESULT bool NextSolution(); +#endif // DEPRECATED: Use TimeLimit() and SetTimeLimit(absl::Duration) instead. // NOTE: These deprecated functions used the convention time_limit = 0 to mean @@ -757,6 +763,7 @@ class MPSolver { friend class SCIPInterface; friend class GurobiInterface; friend class CplexInterface; + friend class XpressInterface; friend class SLMInterface; friend class MPSolverInterface; friend class GLOPInterface; @@ -967,6 +974,7 @@ class MPObjective { friend class SLMInterface; friend class GurobiInterface; friend class CplexInterface; + friend class XpressInterface; friend class GLOPInterface; friend class BopInterface; friend class SatInterface; @@ -1074,6 +1082,7 @@ class MPVariable { friend class SLMInterface; friend class GurobiInterface; friend class CplexInterface; + friend class XpressInterface; friend class GLOPInterface; friend class MPVariableSolutionValueTest; friend class BopInterface; @@ -1215,6 +1224,7 @@ class MPConstraint { friend class SLMInterface; friend class GurobiInterface; friend class CplexInterface; + friend class XpressInterface; friend class GLOPInterface; friend class BopInterface; friend class SatInterface; diff --git a/ortools/linear_solver/xpress_interface.cc b/ortools/linear_solver/xpress_interface.cc new file mode 100644 index 0000000000..90c75f49f3 --- /dev/null +++ b/ortools/linear_solver/xpress_interface.cc @@ -0,0 +1,1485 @@ +// Copyright 2019 RTE +// 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. + +// Initial version of this code was provided by RTE + +#include +#include + +#include "absl/strings/str_format.h" +#include "ortools/base/integral_types.h" +#include "ortools/base/logging.h" +#include "ortools/base/timer.h" +#include "ortools/linear_solver/linear_solver.h" + +#if defined(USE_XPRESS) + +extern "C" { +#include "xprs.h" +} + +#define XPRS_INTEGER 'I' +#define XPRS_CONTINUOUS 'C' + +void printError(const XPRSprob& mLp, int line) { + char errmsg[512]; + XPRSgetlasterror(mLp, errmsg); + VLOG(0) << ("Function line %d did not execute correctly: %s\n", line, errmsg); + exit(0); +} + +int XPRSgetnumcols(const XPRSprob& mLp) +{ + int nCols = 0; + XPRSgetintattrib(mLp, XPRS_COLS, &nCols); + return nCols; +} + +int XPRSgetnumrows(const XPRSprob& mLp) +{ + int nRows = 0; + XPRSgetintattrib(mLp, XPRS_ROWS, &nRows); + return nRows; +} + +int XPRSgetitcnt(const XPRSprob& mLp) +{ + int nIters = 0; + XPRSgetintattrib(mLp, XPRS_SIMPLEXITER, &nIters); + return nIters; +} + +int XPRSgetnodecnt(const XPRSprob& mLp) +{ + int nNodes = 0; + XPRSgetintattrib(mLp, XPRS_NODES, &nNodes); + return nNodes; +} + +int XPRSsetobjoffset(const XPRSprob& mLp, double value) { XPRSsetdblcontrol(mLp, XPRS_OBJRHS, value); return 0; } + +enum XPRS_BASIS_STATUS { + XPRS_AT_LOWER = 0, + XPRS_BASIC = 1, + XPRS_AT_UPPER = 2, + XPRS_FREE_SUPER = 3 +}; + +// In case we need to return a double but don't have a value for that +// we just return a NaN. +#if !defined(CPX_NAN) +#define XPRS_NAN std::numeric_limits::quiet_NaN() +#endif + +// The argument to this macro is the invocation of a XPRS function that +// returns a status. If the function returns non-zero the macro aborts +// the program with an appropriate error message. +#define CHECK_STATUS(s) \ + do { \ + int const status_ = s; \ + CHECK_EQ(0, status_); \ + } while (0) + +namespace operations_research { + + using std::unique_ptr; + + // For a model that is extracted to an instance of this class there is a + // 1:1 corresponence between MPVariable instances and XPRESS columns: the + // index of an extracted variable is the column index in the XPRESS model. + // Similiar for instances of MPConstraint: the index of the constraint in + // the model is the row index in the XPRESS model. + class XpressInterface : public MPSolverInterface { + public: + // NOTE: 'mip' specifies the type of the problem (either continuous or + // mixed integer. This type is fixed for the lifetime of the + // instance. There are no dynamic changes to the model type. + explicit XpressInterface(MPSolver *const solver, bool mip); + ~XpressInterface(); + + // Sets the optimization direction (min/max). + virtual void SetOptimizationDirection(bool maximize); + + // ----- Solve ----- + // Solve the problem using the parameter values specified. + virtual MPSolver::ResultStatus Solve(MPSolverParameters const ¶m); + + // ----- Model modifications and extraction ----- + // Resets extracted model + virtual void Reset(); + + virtual void SetVariableBounds(int var_index, double lb, double ub); + virtual void SetVariableInteger(int var_index, bool integer); + virtual void SetConstraintBounds(int row_index, double lb, double ub); + + virtual void AddRowConstraint(MPConstraint *const ct); + virtual void AddVariable(MPVariable *const var); + virtual void SetCoefficient(MPConstraint *const constraint, + MPVariable const *const variable, + double new_value, double old_value); + + // Clear a constraint from all its terms. + virtual void ClearConstraint(MPConstraint *const constraint); + // Change a coefficient in the linear objective + virtual void SetObjectiveCoefficient(MPVariable const *const variable, + double coefficient); + // Change the constant term in the linear objective. + virtual void SetObjectiveOffset(double value); + // Clear the objective from all its terms. + virtual void ClearObjective(); + + // ------ Query statistics on the solution and the solve ------ + // Number of simplex iterations + virtual int64 iterations() const; + // Number of branch-and-bound nodes. Only available for discrete problems. + virtual int64 nodes() const; + // Best objective bound. Only available for discrete problems. + virtual double best_objective_bound() const; + + // Returns the basis status of a row. + virtual MPSolver::BasisStatus row_status(int constraint_index) const; + // Returns the basis status of a column. + virtual MPSolver::BasisStatus column_status(int variable_index) const; + + // ----- Misc ----- + + // Query problem type. + // Remember that problem type is a static property that is set + // in the constructor and never changed. + virtual bool IsContinuous() const { return IsLP(); } + virtual bool IsLP() const { return !mMip; } + virtual bool IsMIP() const { return mMip; } + + virtual void ExtractNewVariables(); + virtual void ExtractNewConstraints(); + virtual void ExtractObjective(); + + virtual std::string SolverVersion() const; + + virtual void *underlying_solver() { return reinterpret_cast(mLp); } + + virtual double ComputeExactConditionNumber() const { + if (!IsContinuous()) { + LOG(DFATAL) << "ComputeExactConditionNumber not implemented for" + << " XPRESS_MIXED_INTEGER_PROGRAMMING"; + return 0.0; + } + + // TODO(user,user): Not yet working. + LOG(DFATAL) << "ComputeExactConditionNumber not implemented for" + << " XPRESS_LINEAR_PROGRAMMING"; + return 0.0; + } + + protected: + // Set all parameters in the underlying solver. + virtual void SetParameters(MPSolverParameters const ¶m); + // Set each parameter in the underlying solver. + virtual void SetRelativeMipGap(double value); + virtual void SetPrimalTolerance(double value); + virtual void SetDualTolerance(double value); + virtual void SetPresolveMode(int value); + virtual void SetScalingMode(int value); + virtual void SetLpAlgorithm(int value); + + virtual bool ReadParameterFile(std::string const &filename); + virtual std::string ValidFileExtensionForParameterFile() const; + + private: + // Mark modeling object "out of sync". This implicitly invalidates + // solution information as well. It is the counterpart of + // MPSolverInterface::InvalidateSolutionSynchronization + void InvalidateModelSynchronization() { + mCstat = 0; + mRstat = 0; + sync_status_ = MUST_RELOAD; + } + + // Transform XPRESS basis status to MPSolver basis status. + static MPSolver::BasisStatus xformBasisStatus(int xpress_basis_status); + + private: + XPRSprob mLp; + bool const mMip; + // Incremental extraction. + // Without incremental extraction we have to re-extract the model every + // time we perform a solve. Due to the way the Reset() function is + // implemented, this will lose MIP start or basis information from a + // previous solve. On the other hand, if there is a significant changes + // to the model then just re-extracting everything is usually faster than + // keeping the low-level modeling object in sync with the high-level + // variables/constraints. + // Note that incremental extraction is particularly expensive in function + // ExtractNewVariables() since there we must scan _all_ old constraints + // and update them with respect to the new variables. + bool const supportIncrementalExtraction; + + // Use slow and immediate updates or try to do bulk updates. + // For many updates to the model we have the option to either perform + // the update immediately with a potentially slow operation or to + // just mark the low-level modeling object out of sync and re-extract + // the model later. + enum SlowUpdates { + SlowSetCoefficient = 0x0001, + SlowClearConstraint = 0x0002, + SlowSetObjectiveCoefficient = 0x0004, + SlowClearObjective = 0x0008, + SlowSetConstraintBounds = 0x0010, + SlowSetVariableInteger = 0x0020, + SlowSetVariableBounds = 0x0040, + SlowUpdatesAll = 0xffff + } const slowUpdates; + // XPRESS has no method to query the basis status of a single variable. + // Hence we query the status only once and cache the array. This is + // much faster in case the basis status of more than one row/column + // is required. + unique_ptr mutable mCstat; + unique_ptr mutable mRstat; + + // Setup the right-hand side of a constraint from its lower and upper bound. + static void MakeRhs(double lb, double ub, double &rhs, char &sense, + double &range); + }; + + /** init XPRESS environment */ + int init_xpress_env(int xpress_oem_license_key = 0); + int init_xpress_env(int xpress_oem_license_key) + { + int code; + + char *xpresspath; + xpresspath = getenv("XPRESS"); + + google::InitGoogleLogging("Xpress"); + + if (!xpresspath) + { + VLOG(0) << "XpressInterface Error : Environment variable XPRESS undefined.\n"; + return -1; + } + + /** if not an OEM key */ + if (xpress_oem_license_key == 0) + { + VLOG(0) << "XpressInterface : Initialising xpress-MP with parameter " << xpresspath << std::endl; + + code = XPRSinit(xpresspath); + + if (!code) + { + /** XPRSbanner informs about Xpress version, options and error messages */ + char banner[1000]; + XPRSgetbanner(banner); + + VLOG(0) << "XpressInterface : Xpress banner :\n" << banner << std::endl; + + return 0; + } + else + { + char errmsg[256]; + XPRSgetlicerrmsg(errmsg, 256); + + VLOG(0) << "XpressInterface : License error : " << errmsg << std::endl; + VLOG(0) << "XpressInterface : XPRSinit returned code : "<< code <<"\n"; + + char banner[1000]; + XPRSgetbanner(banner); + + VLOG(google::GLOG_ERROR) << "XpressInterface : Xpress banner :\n"<< banner << "\n"; + return -1; + } + } + /** if OEM key */ + else + { + VLOG(0) << "XpressInterface : Initialising xpress-MP with OEM key " << xpress_oem_license_key << "\n"; + + int nvalue = 0; + int ierr; + char slicmsg[256] = ""; + char errmsg[256]; + + XPRSlicense(&nvalue, slicmsg); + VLOG(0) << "XpressInterface : First message from XPRSLicense : " << slicmsg << "\n"; + + nvalue = xpress_oem_license_key - ((nvalue * nvalue) / 19); + ierr = XPRSlicense(&nvalue, slicmsg); + + VLOG(0) << "XpressInterface : Second message from XPRSLicense : " << slicmsg << "\n"; + if (ierr == 16) + { + VLOG(0) << "XpressInterface : Optimizer development software detected\n"; + } + else if (ierr != 0) + { + /** get the license error message */ + XPRSgetlicerrmsg(errmsg, 256); + + VLOG(google::GLOG_ERROR) << "XpressInterface : " << errmsg << "\n"; + return -1; + } + + code = XPRSinit(NULL); + + if (!code) + return 0; + else + { + VLOG(google::GLOG_ERROR) << "XPRSinit returned code : " << code << "\n"; + return -1; + } + } + } + + // Creates a LP/MIP instance. + XpressInterface::XpressInterface(MPSolver *const solver, bool mip) + : MPSolverInterface(solver), + mLp(0), + mMip(mip), + slowUpdates(static_cast(SlowSetObjectiveCoefficient | + SlowClearObjective)), + supportIncrementalExtraction(false), + mCstat(), + mRstat() { + + char const *name = solver_->name_.c_str(); + int status = init_xpress_env(); + CHECK_STATUS(status); + status = XPRScreateprob(&mLp); + CHECK_STATUS(status); + DCHECK(mLp != nullptr); // should not be NULL if status=0 + CHECK_STATUS(XPRSloadlp(mLp, "newProb", 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); + + CHECK_STATUS(XPRSchgobjsense(mLp, maximize_ ? XPRS_OBJ_MAXIMIZE : XPRS_OBJ_MINIMIZE)); + } + + XpressInterface::~XpressInterface() { + CHECK_STATUS(XPRSdestroyprob(mLp)); + google::ShutdownGoogleLogging(); + } + + std::string XpressInterface::SolverVersion() const { + // We prefer XPRSversionnumber() over XPRSversion() since the + // former will never pose any encoding issues. + int version = 0; + CHECK_STATUS(XPRSgetintcontrol(mLp, XPRS_VERSION, &version)); + + int const major = version / 1000000; + version -= major * 1000000; + int const release = version / 10000; + version -= release * 10000; + int const mod = version / 100; + version -= mod * 100; + int const fix = version; + + return absl::StrFormat("XPRESS library version %d.%02d.%02d.%02d", major, release, + mod, fix); + } + + // ------ Model modifications and extraction ----- + + void XpressInterface::Reset() { + // Instead of explicitly clearing all modeling objects we + // just delete the problem object and allocate a new one. + CHECK_STATUS(XPRSdestroyprob(mLp)); + + int status; + const char *const name = solver_->name_.c_str(); + status = XPRScreateprob(&mLp); + CHECK_STATUS(status); + DCHECK(mLp != nullptr); // should not be NULL if status=0 + CHECK_STATUS(XPRSloadlp(mLp, "newProb", 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); + + CHECK_STATUS(XPRSchgobjsense(mLp, maximize_ ? XPRS_OBJ_MAXIMIZE : XPRS_OBJ_MINIMIZE)); + + ResetExtractionInformation(); + mCstat = 0; + mRstat = 0; + } + + void XpressInterface::SetOptimizationDirection(bool maximize) { + InvalidateSolutionSynchronization(); + XPRSchgobjsense(mLp, maximize ? XPRS_OBJ_MAXIMIZE : XPRS_OBJ_MINIMIZE); + } + + void XpressInterface::SetVariableBounds(int var_index, double lb, double ub) { + InvalidateSolutionSynchronization(); + + // Changing the bounds of a variable is fast. However, doing this for + // many variables may still be slow. So we don't perform the update by + // default. However, if we support incremental extraction + // (supportIncrementalExtraction is true) then we MUST perform the + // update here or we will lose it. + + if (!supportIncrementalExtraction && !(slowUpdates & SlowSetVariableBounds)) { + InvalidateModelSynchronization(); + } + else { + if (variable_is_extracted(var_index)) { + // Variable has already been extracted, so we must modify the + // modeling object. + DCHECK_LT(var_index, last_variable_index_); + char const lu[2] = { 'L', 'U' }; + double const bd[2] = { lb, ub }; + int const idx[2] = { var_index, var_index }; + CHECK_STATUS(XPRSchgbounds(mLp, 2, idx, lu, bd)); + } + else { + // Variable is not yet extracted. It is sufficient to just mark + // the modeling object "out of sync" + InvalidateModelSynchronization(); + } + } + } + + // Modifies integrality of an extracted variable. + void XpressInterface::SetVariableInteger(int var_index, bool integer) { + InvalidateSolutionSynchronization(); + + // NOTE: The type of the model (continuous or mixed integer) is + // defined once and for all in the constructor. There are no + // dynamic changes to the model type. + + // Changing the type of a variable should be fast. Still, doing all + // updates in one big chunk right before solve() is usually faster. + // However, if we support incremental extraction + // (supportIncrementalExtraction is true) then we MUST change the + // type of extracted variables here. + + if (!supportIncrementalExtraction && + !(slowUpdates && SlowSetVariableInteger)) { + InvalidateModelSynchronization(); + } + else { + if (mMip) { + if (variable_is_extracted(var_index)) { + // Variable is extracted. Change the type immediately. + // TODO: Should we check the current type and don't do anything + // in case the type does not change? + DCHECK_LE(var_index, XPRSgetnumcols(mLp)); + char const type = integer ? XPRS_INTEGER : XPRS_CONTINUOUS; + CHECK_STATUS(XPRSchgcoltype(mLp, 1, &var_index, &type)); + } + else + InvalidateModelSynchronization(); + } + else { + LOG(DFATAL) + << "Attempt to change variable to integer in non-MIP problem!"; + } + } + } + + // Setup the right-hand side of a constraint. + void XpressInterface::MakeRhs(double lb, double ub, double &rhs, char &sense, + double &range) { + if (lb == ub) { + // Both bounds are equal -> this is an equality constraint + rhs = lb; + range = 0.0; + sense = 'E'; + } + else if (lb > XPRS_MINUSINFINITY && ub < XPRS_PLUSINFINITY) { + // Both bounds are finite -> this is a ranged constraint + // The value of a ranged constraint is allowed to be in + // [ rhs[i], rhs[i]+rngval[i] ] + // see also the reference documentation for XPRSnewrows() + if (ub < lb) { + // The bounds for the constraint are contradictory. XPRESS models + // a range constraint l <= ax <= u as + // ax = l + v + // where v is an auxiliary variable the range of which is controlled + // by l and u: if l < u then v in [0, u-l] + // else v in [u-l, 0] + // (the range is specified as the rngval[] argument to XPRSnewrows). + // Thus XPRESS cannot represent range constraints with contradictory + // bounds and we must error out here. + CHECK_STATUS(-1); + } + rhs = ub; + range = ub - lb; + sense = 'R'; + } + else if (ub < XPRS_PLUSINFINITY || + (std::abs(ub) == XPRS_PLUSINFINITY && std::abs(lb) > XPRS_PLUSINFINITY)) { + // Finite upper, infinite lower bound -> this is a <= constraint + rhs = ub; + range = 0.0; + sense = 'L'; + } + else if (lb > XPRS_MINUSINFINITY || + (std::abs(lb) == XPRS_PLUSINFINITY && std::abs(ub) > XPRS_PLUSINFINITY)) { + // Finite lower, infinite upper bound -> this is a >= constraint + rhs = lb; + range = 0.0; + sense = 'G'; + } + else { + // Lower and upper bound are both infinite. + // This is used for example in .mps files to specify alternate + // objective functions. + // Note that the case lb==ub was already handled above, so we just + // pick the bound with larger magnitude and create a constraint for it. + // Note that we replace the infinite bound by XPRS_PLUSINFINITY since + // bounds with larger magnitude may cause other XPRESS functions to + // fail (for example the export to LP files). + DCHECK_GT(std::abs(lb), XPRS_PLUSINFINITY); + DCHECK_GT(std::abs(ub), XPRS_PLUSINFINITY); + if (std::abs(lb) > std::abs(ub)) { + rhs = (lb < 0) ? -XPRS_PLUSINFINITY : XPRS_PLUSINFINITY; + sense = 'G'; + } + else { + rhs = (ub < 0) ? -XPRS_PLUSINFINITY : XPRS_PLUSINFINITY; + sense = 'L'; + } + range = 0.0; + } + } + + void XpressInterface::SetConstraintBounds(int index, double lb, double ub) { + InvalidateSolutionSynchronization(); + + // Changing rhs, sense, or range of a constraint is not too slow. + // Still, doing all the updates in one large operation is faster. + // Note however that if we do not want to re-extract the full model + // for each solve (supportIncrementalExtraction is true) then we MUST + // update the constraint here, otherwise we lose this update information. + + if (!supportIncrementalExtraction && + !(slowUpdates & SlowSetConstraintBounds)) { + InvalidateModelSynchronization(); + } + else { + if (constraint_is_extracted(index)) { + // Constraint is already extracted, so we must update its bounds + // and its type. + DCHECK(mLp != NULL); + char sense; + double range, rhs; + MakeRhs(lb, ub, rhs, sense, range); + CHECK_STATUS(XPRSchgrhs(mLp, 1, &index, &lb)); + CHECK_STATUS(XPRSchgrowtype(mLp, 1, &index, &sense)); + CHECK_STATUS(XPRSchgrhsrange(mLp, 1, &index, &range)); + } + else { + // Constraint is not yet extracted. It is sufficient to mark the + // modeling object as "out of sync" + InvalidateModelSynchronization(); + } + } + } + + void XpressInterface::AddRowConstraint(MPConstraint *const ct) { + // This is currently only invoked when a new constraint is created, + // see MPSolver::MakeRowConstraint(). + // At this point we only have the lower and upper bounds of the + // constraint. We could immediately call XPRSaddrows() here but it is + // usually much faster to handle the fully populated constraint in + // ExtractNewConstraints() right before the solve. + InvalidateModelSynchronization(); + } + + void XpressInterface::AddVariable(MPVariable *const ct) { + // This is currently only invoked when a new variable is created, + // see MPSolver::MakeVar(). + // At this point the variable does not appear in any constraints or + // the objective function. We could invoke XPRSaddcols() to immediately + // create the variable here but it is usually much faster to handle the + // fully setup variable in ExtractNewVariables() right before the solve. + InvalidateModelSynchronization(); + } + + void XpressInterface::SetCoefficient(MPConstraint *const constraint, + MPVariable const *const variable, + double new_value, double) { + InvalidateSolutionSynchronization(); + + // Changing a single coefficient in the matrix is potentially pretty + // slow since that coefficient has to be found in the sparse matrix + // representation. So by default we don't perform this update immediately + // but instead mark the low-level modeling object "out of sync". + // If we want to support incremental extraction then we MUST perform + // the modification immediately or we will lose it. + + if (!supportIncrementalExtraction && !(slowUpdates & SlowSetCoefficient)) { + InvalidateModelSynchronization(); + } + else { + int const row = constraint->index(); + int const col = variable->index(); + if (constraint_is_extracted(row) && variable_is_extracted(col)) { + // If row and column are both extracted then we can directly + // update the modeling object + DCHECK_LE(row, last_constraint_index_); + DCHECK_LE(col, last_variable_index_); + CHECK_STATUS(XPRSchgcoef(mLp, row, col, new_value)); + } + else { + // If either row or column is not yet extracted then we can + // defer the update to ExtractModel() + InvalidateModelSynchronization(); + } + } + } + + void XpressInterface::ClearConstraint(MPConstraint *const constraint) { + int const row = constraint->index(); + if (!constraint_is_extracted(row)) + // There is nothing to do if the constraint was not even extracted. + return; + + // Clearing a constraint means setting all coefficients in the corresponding + // row to 0 (we cannot just delete the row since that would renumber all + // the constraints/rows after it). + // Modifying coefficients in the matrix is potentially pretty expensive + // since they must be found in the sparse matrix representation. That is + // why by default we do not modify the coefficients here but only mark + // the low-level modeling object "out of sync". + + if (!(slowUpdates & SlowClearConstraint)) { + InvalidateModelSynchronization(); + } + else { + InvalidateSolutionSynchronization(); + + int const len = constraint->coefficients_.size(); + unique_ptr rowind(new int[len]); + unique_ptr colind(new int[len]); + unique_ptr val(new double[len]); + int j = 0; + const auto& coeffs = constraint->coefficients_; + for (auto it(coeffs.begin()); it != coeffs.end(); ++it) { + int const col = it->first->index(); + if (variable_is_extracted(col)) { + rowind[j] = row; + colind[j] = col; + val[j] = 0.0; + ++j; + } + } + if (j) + CHECK_STATUS( + XPRSchgmcoef(mLp, j, rowind.get(), colind.get(), val.get())); + } + } + + void XpressInterface::SetObjectiveCoefficient(MPVariable const *const variable, + double coefficient) { + int const col = variable->index(); + if (!variable_is_extracted(col)) + // Nothing to do if variable was not even extracted + return; + + InvalidateSolutionSynchronization(); + + // The objective function is stored as a dense vector, so updating a + // single coefficient is O(1). So by default we update the low-level + // modeling object here. + // If we support incremental extraction then we have no choice but to + // perform the update immediately. + + if (supportIncrementalExtraction || + (slowUpdates & SlowSetObjectiveCoefficient)) { + CHECK_STATUS(XPRSchgobj(mLp, 1, &col, &coefficient)); + } + else + InvalidateModelSynchronization(); + } + + void XpressInterface::SetObjectiveOffset(double value) { + // Changing the objective offset is O(1), so we always do it immediately. + InvalidateSolutionSynchronization(); + CHECK_STATUS(XPRSsetobjoffset(mLp, value)); + } + + void XpressInterface::ClearObjective() { + InvalidateSolutionSynchronization(); + + // Since the objective function is stored as a dense vector updating + // it is O(n), so we usually perform the update immediately. + // If we want to support incremental extraction then we have no choice + // but to perform the update immediately. + + if (supportIncrementalExtraction || (slowUpdates & SlowClearObjective)) { + int const cols = XPRSgetnumcols(mLp); + unique_ptr ind(new int[cols]); + unique_ptr zero(new double[cols]); + int j = 0; + const auto& coeffs = solver_->objective_->coefficients_; + for (auto it(coeffs.begin()); it != coeffs.end(); ++it) { + int const idx = it->first->index(); + // We only need to reset variables that have been extracted. + if (variable_is_extracted(idx)) { + DCHECK_LT(idx, cols); + ind[j] = idx; + zero[j] = 0.0; + ++j; + } + } + if (j > 0) CHECK_STATUS(XPRSchgobj(mLp, j, ind.get(), zero.get())); + CHECK_STATUS(XPRSsetobjoffset(mLp, 0.0)); + } + else + InvalidateModelSynchronization(); + } + + // ------ Query statistics on the solution and the solve ------ + + int64 XpressInterface::iterations() const { + int iter; + if (!CheckSolutionIsSynchronized()) return kUnknownNumberOfIterations; + return static_cast(XPRSgetitcnt(mLp)); + } + + int64 XpressInterface::nodes() const { + if (mMip) { + if (!CheckSolutionIsSynchronized()) return kUnknownNumberOfNodes; + return static_cast(XPRSgetnodecnt(mLp)); + } + else { + LOG(DFATAL) << "Number of nodes only available for discrete problems"; + return kUnknownNumberOfNodes; + } + } + + // Returns the best objective bound. Only available for discrete problems. + double XpressInterface::best_objective_bound() const { + if (mMip) { + if (!CheckSolutionIsSynchronized() || !CheckBestObjectiveBoundExists()) + // trivial_worst_objective_bound() returns sense*infinity, + // that is meaningful even for infeasible problems + return trivial_worst_objective_bound(); + if (solver_->variables_.size() == 0 && solver_->constraints_.size() == 0) { + // For an empty model the best objective bound is just the offset. + return solver_->Objective().offset(); + } + else { + double value = XPRS_NAN; + CHECK_STATUS(XPRSgetdblattrib(mLp, XPRS_BESTBOUND, &value)); + return value; + } + } + else { + LOG(DFATAL) << "Best objective bound only available for discrete problems"; + return trivial_worst_objective_bound(); + } + } + + // Transform a XPRESS basis status to an MPSolver basis status. + MPSolver::BasisStatus XpressInterface::xformBasisStatus(int xpress_basis_status) { + switch (xpress_basis_status) { + case XPRS_AT_LOWER: + return MPSolver::AT_LOWER_BOUND; + case XPRS_BASIC: + return MPSolver::BASIC; + case XPRS_AT_UPPER: + return MPSolver::AT_UPPER_BOUND; + case XPRS_FREE_SUPER: + return MPSolver::FREE; + default: + LOG(DFATAL) << "Unknown XPRESS basis status"; + return MPSolver::FREE; + } + } + + // Returns the basis status of a row. + MPSolver::BasisStatus XpressInterface::row_status(int constraint_index) const { + if (mMip) { + LOG(FATAL) << "Basis status only available for continuous problems"; + return MPSolver::FREE; + } + + if (CheckSolutionIsSynchronized()) { + if (!mRstat) { + int const rows = XPRSgetnumrows(mLp); + unique_ptr data(new int[rows]); + mRstat.swap(data); + CHECK_STATUS(XPRSgetbasis(mLp, 0, mRstat.get())); + } + } + else + mRstat = 0; + + if (mRstat) + return xformBasisStatus(mRstat[constraint_index]); + else { + LOG(FATAL) << "Row basis status not available"; + return MPSolver::FREE; + } + } + + // Returns the basis status of a column. + MPSolver::BasisStatus XpressInterface::column_status(int variable_index) const { + if (mMip) { + LOG(FATAL) << "Basis status only available for continuous problems"; + return MPSolver::FREE; + } + + if (CheckSolutionIsSynchronized()) { + if (!mCstat) { + int const cols = XPRSgetnumcols(mLp); + unique_ptr data(new int[cols]); + mCstat.swap(data); + CHECK_STATUS(XPRSgetbasis(mLp, mCstat.get(), 0)); + } + } + else + mCstat = 0; + + if (mCstat) + return xformBasisStatus(mCstat[variable_index]); + else { + LOG(FATAL) << "Column basis status not available"; + return MPSolver::FREE; + } + } + + // Extract all variables that have not yet been extracted. + void XpressInterface::ExtractNewVariables() { + // NOTE: The code assumes that a linear expression can never contain + // non-zero duplicates. + + InvalidateSolutionSynchronization(); + + if (!supportIncrementalExtraction) { + // Without incremental extraction ExtractModel() is always called + // to extract the full model. + CHECK(last_variable_index_ == 0 || + last_variable_index_ == solver_->variables_.size()); + CHECK(last_constraint_index_ == 0 || + last_constraint_index_ == solver_->constraints_.size()); + } + + int const last_extracted = last_variable_index_; + int const var_count = solver_->variables_.size(); + int newcols = var_count - last_extracted; + if (newcols > 0) { + // There are non-extracted variables. Extract them now. + + unique_ptr obj(new double[newcols]); + unique_ptr lb(new double[newcols]); + unique_ptr ub(new double[newcols]); + unique_ptr ctype(new char[newcols]); + unique_ptr colname(new const char *[newcols]); + + bool have_names = false; + for (int j = 0, varidx = last_extracted; j < newcols; ++j, ++varidx) { + MPVariable const *const var = solver_->variables_[varidx]; + lb[j] = var->lb(); + ub[j] = var->ub(); + ctype[j] = var->integer() ? XPRS_INTEGER : XPRS_CONTINUOUS; + colname[j] = var->name().empty() ? 0 : var->name().c_str(); + have_names = have_names || var->name().empty(); + obj[j] = solver_->objective_->GetCoefficient(var); + } + + // Arrays for modifying the problem are setup. Update the index + // of variables that will get extracted now. Updating indices + // _before_ the actual extraction makes things much simpler in + // case we support incremental extraction. + // In case of error we just reset the indeces. + std::vector const &variables = solver_->variables(); + for (int j = last_extracted; j < var_count; ++j) { + CHECK(!variable_is_extracted(variables[j]->index())); + set_variable_as_extracted(variables[j]->index(), true); + } + + try { + bool use_newcols = true; + + if (supportIncrementalExtraction) { + // If we support incremental extraction then we must + // update existing constraints with the new variables. + // To do that we use XPRSaddcols() to actually create the + // variables. This is supposed to be faster than combining + // XPRSnewcols() and XPRSchgcoeflist(). + + // For each column count the size of the intersection with + // existing constraints. + unique_ptr collen(new int[newcols]); + for (int j = 0; j < newcols; ++j) collen[j] = 0; + int nonzeros = 0; + // TODO: Use a bitarray to flag the constraints that actually + // intersect new variables? + for (int i = 0; i < last_constraint_index_; ++i) { + MPConstraint const *const ct = solver_->constraints_[i]; + CHECK(constraint_is_extracted(ct->index())); + const auto& coeffs = ct->coefficients_; + for (auto it(coeffs.begin()); it != coeffs.end(); ++it) { + int const idx = it->first->index(); + if (variable_is_extracted(idx) && idx > last_variable_index_) { + collen[idx - last_variable_index_]++; + ++nonzeros; + } + } + } + + if (nonzeros > 0) { + // At least one of the new variables did intersect with an + // old constraint. We have to create the new columns via + // XPRSaddcols(). + use_newcols = false; + unique_ptr begin(new int[newcols + 2]); + unique_ptr cmatind(new int[nonzeros]); + unique_ptr cmatval(new double[nonzeros]); + + // Here is how cmatbeg[] is setup: + // - it is initialized as + // [ 0, 0, collen[0], collen[0]+collen[1], ... ] + // so that cmatbeg[j+1] tells us where in cmatind[] and + // cmatval[] we need to put the next nonzero for column + // j + // - after nonzeros have been setup the array looks like + // [ 0, collen[0], collen[0]+collen[1], ... ] + // so that it is the correct input argument for XPRSaddcols + int *cmatbeg = begin.get(); + cmatbeg[0] = 0; + cmatbeg[1] = 0; + ++cmatbeg; + for (int j = 0; j < newcols; ++j) + cmatbeg[j + 1] = cmatbeg[j] + collen[j]; + + for (int i = 0; i < last_constraint_index_; ++i) { + MPConstraint const *const ct = solver_->constraints_[i]; + int const row = ct->index(); + const auto& coeffs = ct->coefficients_; + for (auto it(coeffs.begin()); it != coeffs.end(); ++it) { + int const idx = it->first->index(); + if (variable_is_extracted(idx) && idx > last_variable_index_) { + int const nz = cmatbeg[idx]++; + cmatind[nz] = row; + cmatval[nz] = it->second; + } + } + } + --cmatbeg; + CHECK_STATUS( + XPRSaddcols(mLp, newcols, nonzeros, obj.get(), cmatbeg, cmatind.get(), cmatval.get(),lb.get(), ub.get()) + ); + } + } + + if (use_newcols) { + // Either incremental extraction is not supported or none of + // the new variables did intersect an existing constraint. + // We can just use XPRSnewcols() to create the new variables. + std::vector collen(newcols, 0); + std::vector cmatbeg(newcols, 0); + unique_ptr cmatind(new int[1]); + unique_ptr cmatval(new double[1]); + cmatind[0] = 0; + cmatval[0] = 1.0; + + CHECK_STATUS( + XPRSaddcols(mLp, newcols, 0, obj.get(), cmatbeg.data(), cmatind.get(), cmatval.get(), lb.get(), ub.get()) + ); + int const cols = XPRSgetnumcols(mLp); + unique_ptr ind(new int[newcols]); + for (int j = 0; j < cols; ++j) + ind[j] = j; + CHECK_STATUS( + XPRSchgcoltype(mLp, cols - last_extracted, ind.get(), ctype.get()) + ); + } + else { + // Incremental extraction: we must update the ctype of the + // newly created variables (XPRSaddcols() does not allow + // specifying the ctype) + if (mMip && XPRSgetnumcols(mLp) > 0) { + // Query the actual number of columns in case we did not + // manage to extract all columns. + int const cols = XPRSgetnumcols(mLp); + unique_ptr ind(new int[newcols]); + for (int j = last_extracted; j < cols; ++j) + ind[j - last_extracted] = j; + CHECK_STATUS(XPRSchgcoltype(mLp, cols - last_extracted, ind.get(), + ctype.get())); + } + } + } + catch (...) { + // Undo all changes in case of error. + int const cols = XPRSgetnumcols(mLp); + if (cols > last_extracted) + { + std::vector colsToDelete; + for (int i = last_extracted; i < cols; ++i) + colsToDelete.push_back(i); + (void)XPRSdelcols(mLp, colsToDelete.size(), colsToDelete.data()); + } + std::vector const &variables = solver_->variables(); + int const size = variables.size(); + for (int j = last_extracted; j < size; ++j) + set_variable_as_extracted(j, false); + throw; + } + } + } + + // Extract constraints that have not yet been extracted. + void XpressInterface::ExtractNewConstraints() { + // NOTE: The code assumes that a linear expression can never contain + // non-zero duplicates. + if (!supportIncrementalExtraction) { + // Without incremental extraction ExtractModel() is always called + // to extract the full model. + CHECK(last_variable_index_ == 0 || + last_variable_index_ == solver_->variables_.size()); + CHECK(last_constraint_index_ == 0 || + last_constraint_index_ == solver_->constraints_.size()); + } + + int const offset = last_constraint_index_; + int const total = solver_->constraints_.size(); + + if (total > offset) { + // There are constraints that are not yet extracted. + + InvalidateSolutionSynchronization(); + + int newCons = total - offset; + int const cols = XPRSgetnumcols(mLp); + DCHECK_EQ(last_variable_index_, cols); + int const chunk = newCons;//10; // max number of rows to add in one shot + + // Update indices of new constraints _before_ actually extracting + // them. In case of error we will just reset the indices. + for (int c = offset; c < total; ++c) + set_constraint_as_extracted(c, true); + + try { + unique_ptr rmatind(new int[cols]); + unique_ptr rmatval(new double[cols]); + unique_ptr rmatbeg(new int[chunk]); + unique_ptr sense(new char[chunk]); + unique_ptr rhs(new double[chunk]); + unique_ptr name(new char const *[chunk]); + unique_ptr rngval(new double[chunk]); + unique_ptr rngind(new int[chunk]); + bool haveRanges = false; + + // Loop over the new constraints, collecting rows for up to + // CHUNK constraints into the arrays so that adding constraints + // is faster. + for (int c = 0; c < newCons; /* nothing */) { + // Collect up to CHUNK constraints into the arrays. + int nextRow = 0; + int nextNz = 0; + for (/* nothing */; c < newCons && nextRow < chunk; ++c, ++nextRow) { + MPConstraint const *const ct = solver_->constraints_[offset + c]; + + // Stop if there is not enough room in the arrays + // to add the current constraint. + if (nextNz + ct->coefficients_.size() > cols) { + DCHECK_GT(nextRow, 0); + break; + } + + // Setup right-hand side of constraint. + MakeRhs(ct->lb(), ct->ub(), rhs[nextRow], sense[nextRow], + rngval[nextRow]); + haveRanges = haveRanges || (rngval[nextRow] != 0.0); + rngind[nextRow] = offset + c; + + // Setup left-hand side of constraint. + rmatbeg[nextRow] = nextNz; + const auto& coeffs = ct->coefficients_; + for (auto it(coeffs.begin()); it != coeffs.end(); ++it) { + int const idx = it->first->index(); + if (variable_is_extracted(idx)) { + DCHECK_LT(nextNz, cols); + DCHECK_LT(idx, cols); + rmatind[nextNz] = idx; + rmatval[nextNz] = it->second; + ++nextNz; + } + } + + // Finally the name of the constraint. + name[nextRow] = ct->name().empty() ? 0 : ct->name().c_str(); + } + if (nextRow > 0) { + CHECK_STATUS( + XPRSaddrows(mLp, nextRow, nextNz, sense.get(), rhs.get(), rngval.get(), rmatbeg.get(), rmatind.get(), rmatval.get()) + ); + if (haveRanges) { + CHECK_STATUS( + XPRSchgrhsrange(mLp, nextRow, rngind.get(), rngval.get())); + } + } + } + } + catch (...) { + // Undo all changes in case of error. + int const rows = XPRSgetnumrows(mLp); + std::vector rowsToDelete; + for (int i = offset; i < rows; ++i) + rowsToDelete.push_back(i); + if (rows > offset) (void)XPRSdelrows(mLp, rowsToDelete.size(), rowsToDelete.data()); + std::vector const &constraints = solver_->constraints(); + int const size = constraints.size(); + for (int i = offset; i < size; ++i) set_constraint_as_extracted(i, false); + throw; + } + } + } + + // Extract the objective function. + void XpressInterface::ExtractObjective() { + // NOTE: The code assumes that the objective expression does not contain + // any non-zero duplicates. + + int const cols = XPRSgetnumcols(mLp); + DCHECK_EQ(last_variable_index_, cols); + + unique_ptr ind(new int[cols]); + unique_ptr val(new double[cols]); + for (int j = 0; j < cols; ++j) { + ind[j] = j; + val[j] = 0.0; + } + + const auto& coeffs = solver_->objective_->coefficients_; + for (auto it = coeffs.begin(); it != coeffs.end(); ++it) { + int const idx = it->first->index(); + if (variable_is_extracted(idx)) { + DCHECK_LT(idx, cols); + val[idx] = it->second; + } + } + + CHECK_STATUS(XPRSchgobj(mLp, cols, ind.get(), val.get())); + CHECK_STATUS(XPRSsetobjoffset(mLp, solver_->Objective().offset())); + } + + // ------ Parameters ----- + + void XpressInterface::SetParameters(const MPSolverParameters ¶m) { + SetCommonParameters(param); + if (mMip) SetMIPParameters(param); + } + + void XpressInterface::SetRelativeMipGap(double value) { + if (mMip) { + CHECK_STATUS(XPRSsetdblcontrol(mLp, XPRS_MIPRELSTOP, value)); + } + else { + LOG(WARNING) << "The relative MIP gap is only available " + << "for discrete problems."; + } + } + + void XpressInterface::SetPrimalTolerance(double value) { + CHECK_STATUS(XPRSsetdblcontrol(mLp, XPRS_FEASTOL, value)); + } + + void XpressInterface::SetDualTolerance(double value) { + CHECK_STATUS(XPRSsetdblcontrol(mLp, XPRS_OPTIMALITYTOL, value)); + } + + void XpressInterface::SetPresolveMode(int value) { + MPSolverParameters::PresolveValues const presolve = + static_cast(value); + + switch (presolve) { + case MPSolverParameters::PRESOLVE_OFF: + CHECK_STATUS(XPRSsetintcontrol(mLp, XPRS_PRESOLVE, 0)); + return; + case MPSolverParameters::PRESOLVE_ON: + CHECK_STATUS(XPRSsetintcontrol(mLp, XPRS_PRESOLVE, 1)); + return; + } + SetIntegerParamToUnsupportedValue(MPSolverParameters::PRESOLVE, value); + } + + // Sets the scaling mode. + void XpressInterface::SetScalingMode(int value) { + MPSolverParameters::ScalingValues const scaling = + static_cast(value); + + switch (scaling) { + case MPSolverParameters::SCALING_OFF: + CHECK_STATUS(XPRSsetintcontrol(mLp, XPRS_SCALING, 0)); + break; + case MPSolverParameters::SCALING_ON: + CHECK_STATUS(XPRSsetintcontrol(mLp, XPRS_SCALING, 1)); + break; + } + } + + // Sets the LP algorithm : primal, dual or barrier. Note that XPRESS offers other + // LP algorithm (e.g. network) and automatic selection + void XpressInterface::SetLpAlgorithm(int value) { + MPSolverParameters::LpAlgorithmValues const algorithm = + static_cast(value); + + int alg = 1; + + switch (algorithm) { + case MPSolverParameters::DUAL: + alg = 2; + break; + case MPSolverParameters::PRIMAL: + alg = 3; + break; + case MPSolverParameters::BARRIER: + alg = 4; + break; + } + + if (alg == XPRS_DEFAULTALG) + SetIntegerParamToUnsupportedValue(MPSolverParameters::LP_ALGORITHM, value); + else { + CHECK_STATUS(XPRSsetintcontrol(mLp, XPRS_DEFAULTALG, alg)); + } + } + + bool XpressInterface::ReadParameterFile(std::string const &filename) { + // Return true on success and false on error. + LOG(DFATAL) << "ReadParameterFile not implemented for XPRESS interface"; + return false; + } + + std::string XpressInterface::ValidFileExtensionForParameterFile() const { + return ".prm"; + } + + MPSolver::ResultStatus XpressInterface::Solve(MPSolverParameters const ¶m) { + int status; + + // Delete chached information + mCstat = 0; + mRstat = 0; + + WallTimer timer; + timer.Start(); + + // Set incrementality + MPSolverParameters::IncrementalityValues const inc = + static_cast( + param.GetIntegerParam(MPSolverParameters::INCREMENTALITY)); + switch (inc) { + case MPSolverParameters::INCREMENTALITY_OFF: + Reset(); /* This should not be required but re-extracting everything + * may be faster, so we do it. */ + break; + case MPSolverParameters::INCREMENTALITY_ON: + XPRSsetintcontrol(mLp, XPRS_CRASH, 0); + break; + } + + // Extract the model to be solved. + // If we don't support incremental extraction and the low-level modeling + // is out of sync then we have to re-extract everything. Note that this + // will lose MIP starts or advanced basis information from a previous + // solve. + if (!supportIncrementalExtraction && sync_status_ == MUST_RELOAD) Reset(); + ExtractModel(); + VLOG(1) << absl::StrFormat("Model build in %.3f seconds.", timer.Get()); + + // Set log level. + CHECK_STATUS( + XPRSsetintcontrol(mLp, XPRS_LPLOG, quiet() ? 0 : 1)); + CHECK_STATUS( + XPRSsetintcontrol(mLp, XPRS_MIPLOG, quiet() ? 0 : 1)); + + // Set parameters. + // NOTE: We must invoke SetSolverSpecificParametersAsString() _first_. + // Its current implementation invokes ReadParameterFile() which in + // turn invokes XPRSreadcopyparam(). The latter will _overwrite_ + // all current parameter settings in the environment. + solver_->SetSolverSpecificParametersAsString( + solver_->solver_specific_parameter_string_); + SetParameters(param); + if (solver_->time_limit()) { + VLOG(1) << "Setting time limit = " << solver_->time_limit() << " ms."; + CHECK_STATUS( + XPRSsetdblcontrol(mLp, XPRS_MAXTIME, solver_->time_limit() * 1e-3)); + } + + // Solve. + // Do not CHECK_STATUS here since some errors (for example CPXERR_NO_MEMORY) + // still allow us to query useful information. + timer.Restart(); + + int xpressstat = 0; + if (mMip) { + if (this->maximize_) + status = XPRSmaxim(mLp, "g"); + else + status = XPRSminim(mLp, "g"); + XPRSgetintattrib(mLp, XPRS_MIPSTATUS, &xpressstat); + } + else { + if (this->maximize_) + status = XPRSmaxim(mLp, ""); + else + status = XPRSminim(mLp, ""); + XPRSgetintattrib(mLp, XPRS_LPSTATUS, &xpressstat); + } + + // Disable screen output right after solve + XPRSsetintcontrol(mLp, XPRS_LPLOG, 0); + + if (status) { + VLOG(1) << absl::StrFormat("Failed to optimize MIP. Error %d", status); + // NOTE: We do not return immediately since there may be information + // to grab (for example an incumbent) + } + else { + VLOG(1) << absl::StrFormat("Solved in %.3f seconds.", timer.Get()); + } + + VLOG(1) << absl::StrFormat("XPRESS solution status %d.", xpressstat); + + // Figure out what solution we have. + bool const feasible = !(xpressstat == XPRS_MIP_INFEAS || xpressstat == XPRS_LP_INFEAS); + + // Get problem dimensions for solution queries below. + int const rows = XPRSgetnumrows(mLp); + int const cols = XPRSgetnumcols(mLp); + DCHECK_EQ(rows, solver_->constraints_.size()); + DCHECK_EQ(cols, solver_->variables_.size()); + + // Capture objective function value. + objective_value_ = XPRS_NAN; + if (feasible) { + if (mMip) { + CHECK_STATUS(XPRSgetdblattrib(mLp, XPRS_MIPOBJVAL, &objective_value_)); + } + else { + CHECK_STATUS(XPRSgetdblattrib(mLp, XPRS_LPOBJVAL, &objective_value_)); + } + } + VLOG(1) << "objective = " << objective_value_; + + // Capture primal and dual solutions + if (mMip) { + // If there is a primal feasible solution then capture it. + if (feasible) { + if (cols > 0) { + unique_ptr x(new double[cols]); + CHECK_STATUS(XPRSgetmipsol(mLp, x.get(), 0)); + for (int i = 0; i < solver_->variables_.size(); ++i) { + MPVariable *const var = solver_->variables_[i]; + var->set_solution_value(x[i]); + VLOG(3) << var->name() << ": value =" << x[i]; + } + } + } + else { + for (int i = 0; i < solver_->variables_.size(); ++i) + solver_->variables_[i]->set_solution_value(XPRS_NAN); + } + + // MIP does not have duals + for (int i = 0; i < solver_->variables_.size(); ++i) + solver_->variables_[i]->set_reduced_cost(XPRS_NAN); + for (int i = 0; i < solver_->constraints_.size(); ++i) + solver_->constraints_[i]->set_dual_value(XPRS_NAN); + } + else { + // Continuous problem. + if (cols > 0) { + unique_ptr x(new double[cols]); + unique_ptr dj(new double[cols]); + if (feasible) CHECK_STATUS(XPRSgetlpsol(mLp, x.get(), 0, 0, dj.get())); + for (int i = 0; i < solver_->variables_.size(); ++i) { + MPVariable *const var = solver_->variables_[i]; + var->set_solution_value(x[i]); + bool value = false, dual = false; + + if (feasible) { + var->set_solution_value(x[i]); + value = true; + } + else + var->set_solution_value(XPRS_NAN); + if (feasible) { + var->set_reduced_cost(dj[i]); + dual = true; + } + else + var->set_reduced_cost(XPRS_NAN); + VLOG(3) << var->name() << ":" + << (value ? absl::StrFormat(" value = %f", x[i]) : "") + << (dual ? absl::StrFormat(" reduced cost = %f", dj[i]) : ""); + } + } + + if (rows > 0) { + unique_ptr pi(new double[rows]); + if (feasible) { + CHECK_STATUS(XPRSgetlpsol(mLp, 0, 0, pi.get(), 0)); + } + for (int i = 0; i < solver_->constraints_.size(); ++i) { + MPConstraint *const ct = solver_->constraints_[i]; + bool dual = false; + if (feasible) { + ct->set_dual_value(pi[i]); + dual = true; + } + else + ct->set_dual_value(XPRS_NAN); + VLOG(4) << "row " << ct->index() << ":" + << (dual ? absl::StrFormat(" dual = %f", pi[i]) : ""); + } + } + } + + // Map XPRESS status to more generic solution status in MPSolver + if (mMip) { + switch (xpressstat) { + case XPRS_MIP_OPTIMAL: + result_status_ = MPSolver::OPTIMAL; + break; + case XPRS_MIP_INFEAS: + result_status_ = MPSolver::INFEASIBLE; + break; + case XPRS_MIP_UNBOUNDED: + result_status_ = MPSolver::UNBOUNDED; + break; + default: + result_status_ = feasible ? MPSolver::FEASIBLE : MPSolver::ABNORMAL; + break; + } + } + else { + switch (xpressstat) { + case XPRS_LP_OPTIMAL: + result_status_ = MPSolver::OPTIMAL; + break; + case XPRS_LP_INFEAS: + result_status_ = MPSolver::INFEASIBLE; + break; + case XPRS_LP_UNBOUNDED: + result_status_ = MPSolver::UNBOUNDED; + break; + default: + result_status_ = feasible ? MPSolver::FEASIBLE : MPSolver::ABNORMAL; + break; + } + } + + sync_status_ = SOLUTION_SYNCHRONIZED; + return result_status_; + } + + MPSolverInterface *BuildXpressInterface(bool mip, MPSolver *const solver) { + return new XpressInterface(solver, mip); + } + +} // namespace operations_research +#endif // #if defined(USE_XPRESS)