diff --git a/solvers/CMakeLists.txt b/solvers/CMakeLists.txt index ad8cfe979..8510de857 100644 --- a/solvers/CMakeLists.txt +++ b/solvers/CMakeLists.txt @@ -103,7 +103,7 @@ function(add_ampl_solver name) if (TARGET ${module}-library) get_target_property(TRG_PROP_IMPRT_LOC ${module}-library IMPORTED_LOCATION) - if (NOT libraries OR TRG_PROP_IMPRT_LOC) + if (NOT libraries)#OR TRG_PROP_IMPRT_LOC) set(libraries ${module}-library) endif() endif () @@ -245,6 +245,9 @@ add_ampl_backend(cbcmp DLL_RUNTIME SHARED_LIB MODULE CBC LIBRARIES ${CBC_LIBS} ${CMAKE_DL_LIBS}) +add_ampl_backend(scip DLL_RUNTIME SHARED_LIB MODULE SCIP + LIBRARIES ${SCIP_LIBS} ${CMAKE_DL_LIBS}) + add_ampl_backend(highsmp DLL_RUNTIME SHARED_LIB MODULE HIGHS LIBRARIES ${HIGHS_LIBS} ${CMAKE_DL_LIBS}) diff --git a/solvers/scip/CHANGES.scip.md b/solvers/scip/CHANGES.scip.md new file mode 100644 index 000000000..f0e76434b --- /dev/null +++ b/solvers/scip/CHANGES.scip.md @@ -0,0 +1,5 @@ +Summary of recent updates to SCIP for AMPL +========================================== + +## Unreleased +- First release of mock driver diff --git a/solvers/scip/README.scip.txt b/solvers/scip/README.scip.txt new file mode 100644 index 000000000..6907aedbc --- /dev/null +++ b/solvers/scip/README.scip.txt @@ -0,0 +1,34 @@ +SCIP driver for AMPL +==================== + +SCIP (Solving Constraint Integer Programs) is currently one of the +fastest non-commercial solvers for mixed integer programming (MIP) +and mixed integer nonlinear programming (MINLP). It is also a framework +for constraint integer programming and branch-cut-and-price. +More information can be found in https://scipopt.org/#scipoptsuite. + +Normally SCIP is invoked by AMPL's solve command, which gives the +invocation + + scip stub -AMPL + +in which stub.nl is an AMPL generic output file (possibly written +by "ampl -obstub" or "ampl -ogstub"). After solving the problem, +scip writes a stub.sol file for use by AMPL's solve and solution +commands. When you run ampl, this all happens automatically if you +give the AMPL commands + + option solver scip; + solve; + +You can control scip either by setting the environment variable +scip_options appropriately (either by using ampl's option command, +or by using the shell's set and export commands before you invoke ampl), +or by passing the options on the command line: + + scip stub [-AMPL] option1=value option2=value ... + +You can put one or more (white-space separated) phrases in $scip_options. +To see the possibilities, invoke + + scip -= diff --git a/solvers/scip/main.cc b/solvers/scip/main.cc new file mode 100644 index 000000000..b558be280 --- /dev/null +++ b/solvers/scip/main.cc @@ -0,0 +1,9 @@ +#include "mp/backend-app.h" + +/// Declare a backend factory +std::unique_ptr CreateScipBackend(); + +extern "C" int main1(int, char **argv) { + return + mp::RunBackendApp(argv, CreateScipBackend); +} diff --git a/solvers/scip/model-mgr-with-std-pb.cc b/solvers/scip/model-mgr-with-std-pb.cc new file mode 100644 index 000000000..2ffcce081 --- /dev/null +++ b/solvers/scip/model-mgr-with-std-pb.cc @@ -0,0 +1,8 @@ +/** + * Generate ModelManagerWithPB + * + * Having a separate .cc should improve compilation speed + */ + +#include "mp/model-mgr-with-std-pb.hpp" + diff --git a/solvers/scip/scip-ampls-c-api.h b/solvers/scip/scip-ampls-c-api.h new file mode 100644 index 000000000..127620fd6 --- /dev/null +++ b/solvers/scip/scip-ampls-c-api.h @@ -0,0 +1,33 @@ +#ifndef SCIPAMPLSCAPI_H +#define SCIPAMPLSCAPI_H +/* + * C API for MP/Scip + */ + +//#include "scip.h" + +#include "mp/ampls-c-api.h" + +/* + * Below are Scip-specific AMPLS API functions. + * They complement the 'public' AMPLS API defined in ampls-c-api.h. + */ + +/// Initialize AMPLS scip. + +/// @param slv_opt: a string of solver options +/// (normally provided in the _options string). +/// Can be NULL. +/// @return pointer to struct AMPLS_MP_Solver to be populated. +void* AMPLSOpenScip(const char* slv_opt, CCallbacks cb); + +/// Shut down solver instance +void AMPLSCloseScip(AMPLS_MP_Solver* slv); + +/// Extract the Scip model handle +void* GetScipmodel(AMPLS_MP_Solver* slv); + + +#endif // SCIPAMPLSCAPI_H + + diff --git a/solvers/scip/scip-lib.c b/solvers/scip/scip-lib.c new file mode 100644 index 000000000..20736312f --- /dev/null +++ b/solvers/scip/scip-lib.c @@ -0,0 +1,30 @@ +#include "scip/scip-ampls-c-api.h" + +#ifdef _WIN32 +#define APIEXPORT __declspec(dllexport) +#else +#define APIEXPORT __attribute__((visibility("default"))) +#endif + + +APIEXPORT void* AMPLloadmodel(int argc, char** argv, CCallbacks cb) { + const char* nl_filename = argv[1]; + const char *slv_opt= argv[2]; + AMPLS_MP_Solver* slv; + slv = AMPLSOpenScip(slv_opt, cb); + if (!slv) + return NULL; + AMPLSLoadNLModel(slv, nl_filename); + return slv; +} +APIEXPORT void* AMPLgetScipmodel(void* slv) { + return GetScipmodel(slv); +} + +APIEXPORT void AMPLwritesolution(AMPLS_MP_Solver* slv, const char* solFileName) { + AMPLSReportResults(slv, solFileName); +} + +APIEXPORT void AMPLclosesolver(AMPLS_MP_Solver* slv) { + AMPLSCloseScip(slv); +} \ No newline at end of file diff --git a/solvers/scip/scip-modelapi-connect.cc b/solvers/scip/scip-modelapi-connect.cc new file mode 100644 index 000000000..77641e7a5 --- /dev/null +++ b/solvers/scip/scip-modelapi-connect.cc @@ -0,0 +1,18 @@ +#include "mp/flat/redef/MIP/converter_mip.h" +#include "mp/flat/model_api_connect.h" + +#include "scipmodelapi.h" + + +namespace mp { + +/// Defining the function in ...-modelapi-connect.cc +/// for recompilation speed +std::unique_ptr +CreateScipModelMgr(ScipCommon& cc, Env& e, + pre::BasicValuePresolver*& pPre) { + return CreateModelMgrWithFlatConverter< + ScipModelAPI, MIPFlatConverter >(cc, e, pPre); +} + +} // namespace mp diff --git a/solvers/scip/scipbackend.cc b/solvers/scip/scipbackend.cc new file mode 100644 index 000000000..3f0315340 --- /dev/null +++ b/solvers/scip/scipbackend.cc @@ -0,0 +1,933 @@ +#include +#include +#include + +#include "mp/env.h" +#include "mp/flat/model_api_base.h" +#include "scipbackend.h" + +extern "C" { + #include "scip-ampls-c-api.h" // Scip AMPLS C API +} +#include "mp/ampls-cpp-api.h" + +namespace { + + +bool InterruptScip(void* scip) { + SCIPinterruptSolve(static_cast(scip)); + return true; +} + +} // namespace {} + +std::unique_ptr CreateScipBackend() { + return std::unique_ptr{new mp::ScipBackend()}; +} + + +namespace mp { + +/// Create Scip Model Manager +/// @param gc: the Scip common handle +/// @param e: environment +/// @param pre: presolver to be returned, +/// need it to convert solution data +/// @return ScipModelMgr +std::unique_ptr +CreateScipModelMgr(ScipCommon&, Env&, pre::BasicValuePresolver*&); + + +ScipBackend::ScipBackend() { + OpenSolver(); + + /// Create a ModelManager + pre::BasicValuePresolver* pPre; + auto data = CreateScipModelMgr(*this, *this, pPre); + SetMM( std::move( data ) ); + SetValuePresolver(pPre); + + /// Copy env/lp to ModelAPI + copy_common_info_to_other(); +} + +ScipBackend::~ScipBackend() { + CloseSolver(); +} + + +const char* ScipBackend::GetBackendName() + { return "SCIPBackend"; } + +std::string ScipBackend::GetSolverVersion() { + return fmt::format("{}.{}.{}", SCIPmajorVersion(), SCIPminorVersion(), + SCIPtechVersion()); +} + + +bool ScipBackend::IsMIP() const { + return SCIPgetNOrigIntVars(getSCIP()) > 0 || SCIPgetNOrigBinVars(getSCIP()) > 0; +} + +bool ScipBackend::IsQCP() const { + for (int i = 0; i < SCIPgetNOrigConss(getSCIP()); i++) { + SCIP_Bool isquadratic; + SCIP_CCALL( SCIPcheckQuadraticNonlinear(getSCIP(), SCIPgetOrigConss(getSCIP())[i], &isquadratic) ); + if (isquadratic == true) + return true; + } + return false; +} + +ArrayRef ScipBackend::PrimalSolution() { + SCIP* scip = getSCIP(); + int num_vars = NumVars(); + std::vector x(num_vars); + for (int i = 0; i < num_vars; i++) + x[i] = SCIPgetSolVal(scip, SCIPgetBestSol(scip), getPROBDATA()->vars[i]); + return x; +} + +pre::ValueMapDbl ScipBackend::DualSolution() { + return {{ { CG_Linear, DualSolution_LP() } }}; +} + +ArrayRef ScipBackend::DualSolution_LP() { + int num_cons = NumLinCons(); + std::vector pi(num_cons); + for (int i=0; i < num_cons; i++) + pi[i] = SCIPgetDualsolLinear(getSCIP(), getPROBDATA()->linconss[i]); + return pi; +} + +double ScipBackend::ObjectiveValue() const { + return SCIPgetPrimalbound(getSCIP()); +} + +double ScipBackend::NodeCount() const { + return SCIPgetNNodes(getSCIP()); +} + +double ScipBackend::SimplexIterations() const { + return SCIPgetNPrimalLPIterations(getSCIP()) + SCIPgetNDualLPIterations(getSCIP()); +} + +int ScipBackend::BarrierIterations() const { + return SCIPgetNBarrierLPIterations(getSCIP()); +} + +void ScipBackend::ExportModel(const std::string &file) { + SCIP_CCALL( SCIPwriteOrigProblem(getSCIP(), file.data(), NULL, FALSE) ); +} + + +void ScipBackend::SetInterrupter(mp::Interrupter *inter) { + inter->SetHandler(InterruptScip, getSCIP()); +} + +void ScipBackend::Solve() { + if (!storedOptions_.exportFile_.empty()) + ExportModel(storedOptions_.exportFile_); + if (!storedOptions_.paramRead_.empty()) + SCIP_CCALL( SCIPreadParams(getSCIP(), storedOptions_.paramRead_.c_str()) ); + if (storedOptions_.heuristics_ != SCIP_PARAMSETTING_DEFAULT) + SCIP_CCALL( SCIPsetHeuristics(getSCIP(), (SCIP_PARAMSETTING)storedOptions_.heuristics_, FALSE) ); + if (storedOptions_.cuts_ != SCIP_PARAMSETTING_DEFAULT) + SCIP_CCALL( SCIPsetSeparating(getSCIP(), (SCIP_PARAMSETTING)storedOptions_.cuts_, FALSE) ); + if (storedOptions_.presolvings_ != SCIP_PARAMSETTING_DEFAULT) + SCIP_CCALL( SCIPsetSeparating(getSCIP(), (SCIP_PARAMSETTING)storedOptions_.presolvings_, FALSE) ); + + + if (storedOptions_.concurrent_) + SCIP_CCALL( SCIPsolveConcurrent(getSCIP()) ); + else + SCIP_CCALL( SCIPsolve(getSCIP()) ); + + WindupSCIPSolve(); +} + +void ScipBackend::WindupSCIPSolve() { } + +void ScipBackend::ReportResults() { + ReportSCIPResults(); + BaseBackend::ReportResults(); +} + +void ScipBackend::ReportSCIPResults() { + SetStatus( ConvertSCIPStatus() ); + AddSCIPMessages(); + if (need_multiple_solutions()) + ReportSCIPPool(); +} +std::vector ScipBackend::getPoolSolution(int i) +{ + SCIP* scip = getSCIP(); + int num_vars = NumVars(); + std::vector vars(num_vars); + for (int j = 0; j < num_vars; j++) + vars[j] = SCIPgetSolVal(scip, SCIPgetSols(scip)[i], getPROBDATA()->vars[j]); + return vars; +} +double ScipBackend::getPoolObjective(int i) +{ + assert(i < SCIPgetNSols(getSCIP())); + double obj; + obj = SCIPgetSolOrigObj(getSCIP(), SCIPgetSols(getSCIP())[i]); + return obj; +} +void ScipBackend::ReportSCIPPool() { + if (!IsMIP()) + return; + int iPoolSolution = -1; + int nsolutions = SCIPgetNSols(getSCIP()); + + while (++iPoolSolution < nsolutions) { + ReportIntermediateSolution( + { getPoolSolution(iPoolSolution), + {}, { getPoolObjective(iPoolSolution) } }); + } +} + + +void ScipBackend::AddSCIPMessages() { + AddToSolverMessage( + fmt::format("{} simplex iterations\n", SimplexIterations())); + if (auto nbi = BarrierIterations()) + AddToSolverMessage( + fmt::format("{} barrier iterations\n", nbi)); + if (!IsContinuous()) { + auto nnd = NodeCount(); + AddToSolverMessage( + fmt::format("{} branching nodes\n", nnd)); + } +} + +std::pair ScipBackend::ConvertSCIPStatus() { + namespace sol = mp::sol; + SCIP_STATUS status = SCIPgetStatus(getSCIP()); + switch (status) { + case SCIP_STATUS_UNKNOWN: + return { sol::UNKNOWN, "solving status not yet known" }; + case SCIP_STATUS_USERINTERRUPT: + return { sol::INTERRUPTED, "unfinished" }; + case SCIP_STATUS_NODELIMIT: + return { sol::LIMIT, "node limit reached" }; + case SCIP_STATUS_TOTALNODELIMIT: + return { sol::LIMIT, "total node limit reached (incl. restarts)" }; + case SCIP_STATUS_STALLNODELIMIT: + return { sol::LIMIT, "stalling node limit reached (no inprovement w.r.t. primal bound)" }; + case SCIP_STATUS_TIMELIMIT: + return { sol::LIMIT, "time limit reached" }; + case SCIP_STATUS_MEMLIMIT: + return { sol::LIMIT, "memory limit reached" }; + case SCIP_STATUS_GAPLIMIT: + return { sol::LIMIT, "gap limit reached" }; + case SCIP_STATUS_SOLLIMIT: + return { sol::LIMIT, "solution limit reached" }; + case SCIP_STATUS_BESTSOLLIMIT: + return { sol::LIMIT, "solution improvement limit reached" }; + case SCIP_STATUS_RESTARTLIMIT: + return { sol::LIMIT, "restart limit was reached" }; + case SCIP_STATUS_OPTIMAL: + return { sol::SOLVED, "optimal solution" }; + case SCIP_STATUS_INFEASIBLE: + return { sol::INFEASIBLE, "infeasible problem" }; + case SCIP_STATUS_UNBOUNDED: + return { sol::UNBOUNDED, "unbounded problem" }; + case SCIP_STATUS_INFORUNBD: + return { sol::INF_OR_UNB, "infeasible or unbounded problem" }; + case SCIP_STATUS_TERMINATE: + return { sol::FAILURE, "process received a SIGTERM signal" }; + } + return { sol::UNKNOWN, "not solved" }; +} + + +void ScipBackend::FinishOptionParsing() { + int v=-1; + GetSolverOption("display/verblevel", v); + set_verbose_mode(v>0); +} + + +////////////////////////////// OPTIONS ///////////////////////////////// +static const mp::OptionValueInfo lp_values_method[] = { + { "s", "automatic simplex (default)", -1}, + { "p", "primal simplex", 0}, + { "d", "dual simplex", 1}, + { "b", "barrier", 2}, + { "c", "barrier with crossover", 3}, +}; + +static const mp::OptionValueInfo values_pricing[] = { + { "l", "lpi default (default)", -1}, + { "a", "auto", 0}, + { "f", "full pricing", 1}, + { "p", "partial", 2}, + { "s", "steepest edge pricing", 3}, + { "q", "quickstart steepest edge pricing", 4}, + { "d", "devex pricing", 5} +}; + +static const mp::OptionValueInfo estimation_method[] = { + { "c", "completion", 0}, + { "e", "ensemble", 1}, + { "g", "time series forecasts on either gap", 2}, + { "l", "leaf frequency", 3}, + { "o", "open nodes", 4}, + { "w", "tree weight (default)", 5}, + { "s", "ssg", 6}, + { "t", "tree profile", 7}, + { "b", "wbe ", 8} +}; + +static const mp::OptionValueInfo estimation_completion[] = { + { "a", "auto (default)", 0}, + { "g", "gap", 1}, + { "w", "tree weight", 2}, + { "m", "monotone regression", 3}, + { "r", "regression forest", 4}, + { "s", "ssg", 5} +}; + +static const mp::OptionValueInfo childsel[] = { + { "d", "down", 0}, + { "u", "up", 1}, + { "p", "pseudo costs", 2}, + { "i", "inference", 3}, + { "l", "lp value", 4}, + { "r", "root LP value difference", 5}, + { "h", "hybrid inference/root LP value difference (default)", 6} +}; + + +void ScipBackend::InitCustomOptions() { + + set_option_header( + "SCIP Optimizer Options for AMPL\n" + "--------------------------------------------\n" + "\n" + "To set these options, assign a string specifying their values to the " + "AMPL option ``scip_options``. For example::\n" + "\n" + " ampl: option scip_options 'mipgap=1e-6';\n"); + + AddSolverOption("tech:outlev outlev", + "0*/1/2/3/4/5: Whether to write SCIP log lines (chatter) to stdout and to file.", + "display/verblevel", 0, 5); + + AddStoredOption("tech:exportfile writeprob writemodel", + "Specifies the name of a file where to export the model before " + "solving it. This file name can have extension ``.lp``, ``.mps``, etc. " + "Default = \"\" (don't export the model).", + storedOptions_.exportFile_); + + AddStoredOption("tech:logfile logfile", + "Log file name; note that the solver log will be written to the log " + "regardless of the value of tech:outlev.", + storedOptions_.logFile_); + + AddStoredOption("tech:param:read param:read paramfile", + "Filename of SCIP parameter file (as path)." + "The suffix on a parameter file should be .set.\n", + storedOptions_.paramRead_); + + AddSolverOption("alg:method method lpmethod", + "LP algorithm for solving initial LP relaxations:\n" + "\n.. value-table::\n", "lp/initalgorithm", lp_values_method, "s"); + + AddSolverOption("alg:remethod remethod relpmethod", + "LP algorithm for resolving LP relaxations if a starting basis exists:\n" + "\n.. value-table::\n", "lp/resolvealgorithm", lp_values_method, "s"); + + AddStoredOption("alg:concurrent concurrent", + "0/1: whether to solve the problem using concurrent solvers" + "\n" + " | 0 - No concurrent solvers are used to solve the problem (default)\n" + " | 1 - Concurrent solvers are used to solve the problem.", + storedOptions_.concurrent_, 0, 1); + + /////////////////////// BRANCHING /////////////////////// + AddSolverOption("branch:preferbinary preferbinary", + "0/1: whether branching on binary variables should be preferred" + "\n" + " | 0 - Binary variables should not be preferred on branching (default)\n" + " | 1 - Binary variables should be preferred on branching.", + "branching/preferbinary", 0, 1); + + //////////////////////// CUTS ////////////////////////// + AddSolverOption("cut:dircutoffdistweight dircutoffdistweight", + "Weight of directed cutoff distance in cut score calculation (default: 0.0)", + "cutselection/hybrid/dircutoffdistweight", 0.0, SCIP_REAL_MAX); + + AddSolverOption("cut:efficacyweight efficacyweight", + "Weight of efficacy in cut score calculation (default: 1.0)", + "cutselection/hybrid/efficacyweight", 0.0, SCIP_REAL_MAX); + + AddSolverOption("cut:intsupportweight intsupportweight", + "Weight of integral support in cut score calculation (default: 0.1)", + "cutselection/hybrid/intsupportweight", 0.0, SCIP_REAL_MAX); + + AddSolverOption("cut:minortho minortho", + "Minimal orthogonality for a cut to enter the LP (default: 0.9)", + "cutselection/hybrid/minortho", 0.0, SCIP_REAL_MAX); + + AddSolverOption("cut:minorthoroot minorthoroot", + "Minimal orthogonality for a cut to enter the LP in the root node (default: 0.1)", + "cutselection/hybrid/minorthoroot", 0.0, SCIP_REAL_MAX); + + AddSolverOption("cut:objparalweight objparalweight", + "Weight of objective parallelism in cut score calculation (default: 0.1)", + "cutselection/hybrid/objparalweight", 0.0, SCIP_REAL_MAX); + + + AddSolverOption("cut:maxcuts maxcuts", + "Maximal number of cuts separated per separation round (0: disable local separation; default: 100)", + "separating/maxcuts", 0, INT_MAX); + + AddSolverOption("cut:maxcutsroot maxcutsroot", + "Maximal number of separated cuts at the root node (0: disable root node separation; default: 5000)", + "separating/maxcutsroot", 0, INT_MAX); + + AddSolverOption("cut:maxrounds", + "Maximal number of separation rounds per node (default: -1: unlimited)", + "separating/maxrounds", -1, INT_MAX); + + AddSolverOption("cut:maxroundsroot", + "Maximal number of separation rounds in the root node (default: -1: unlimited)", + "separating/maxroundsroot", -1, INT_MAX); + + AddSolverOption("cut:maxstallrounds", + "Maximal number of consecutive separation rounds without objective or integrality improvement in local nodes (-1: no additional restriction; default: 1)", + "separating/maxstallrounds", -1, INT_MAX); + + AddSolverOption("cut:maxstallroundsroot", + "Maximal number of consecutive separation rounds without objective or integrality improvement in the root node (-1: no additional restriction; default: 10)", + "separating/maxstallroundsroot", -1, INT_MAX); + + AddSolverOption("cut:minactivityquot", + "Minimum cut activity quotient to convert cuts into constraints during a restart (0.0: all cuts are converted; default: 0.8)", + "separating/minactivityquot", 0.0, 1.0); + + AddSolverOption("cut:minefficacy", + "Minimal efficacy for a cut to enter the LP (default: 0.0001)", + "separating/minefficacy", 0.0, SCIP_REAL_MAX); + + AddSolverOption("cut:minefficacyroot", + "Minimal efficacy for a cut to enter the LP in the root node (default: 0.0001)", + "separating/minefficacyroot", 0.0, SCIP_REAL_MAX); + + AddSolverOption("cut:poolfreq", + "Separation frequency for the global cut pool (-1: disable global cut pool; 0: only separate pool at the root; default: 10)", + "separating/poolfreq", -1, SCIP_MAXTREEDEPTH); + + AddStoredOption("cut:settings", + "0/1/2/3: sets cuts settings" + "\n" + " | 0 - Sets cuts default (default)\n" + " | 1 - Sets cuts aggressive\n" + " | 2 - Sets cuts fast\n" + " | 3 - Sets cuts off.", + storedOptions_.cuts_, 0, 3); + + ////////////////////// HEURISTICS ////////////////////// + AddStoredOption("heu:settings", + "0/1/2/3: sets heuristics settings" + "\n" + " | 0 - Sets heuristics default (default)\n" + " | 1 - Sets heuristics aggressive\n" + " | 2 - Sets heuristics fast\n" + " | 3 - Sets heuristics off.", + storedOptions_.heuristics_, 0, 3); + + //////////////////////// LIMITS //////////////////////// + AddSolverOption("lim:absgap absgap", + "Solving stops, if the absolute gap = |primalbound - dualbound| is below the given value (default: 0.0)", + "limits/absgap", 0.0, SCIP_REAL_MAX); + + AddSolverOption("lim:autorestartnodes", + "If solve exceeds this number of nodes for the first time, an automatic restart is triggered (default: -1: no automatic restart)", + "limits/autorestartnodes", -1, INT_MAX); + + AddSolverOption("lim:bestsol", + "Solving stops, if the given number of solution improvements were found (default: -1: no limit)", + "limits/bestsol", -1, INT_MAX); + + AddSolverOption("lim:gap gap", + "Solving stops, if the relative gap = |primal - dual|/MIN(|dual|,|primal|) is below the given value, the gap is 'Infinity', if primal and dual bound have opposite signs (default: 0.0)", + "limits/gap", 0.0, SCIP_REAL_MAX); + + AddSolverOption("lim:maxorigsol", + "Maximal number of solutions candidates to store in the solution storage of the original problem (default: 10)", + "limits/maxorigsol", 0, INT_MAX); + + AddSolverOption("lim:maxsol", + "Maximal number of solutions to store in the solution storage (default: 100)", + "limits/maxsol", 1, INT_MAX); + + AddSolverOption("lim:memory memory", + "#maximal memory usage in MB; reported memory usage is lower than real memory usage! (default: 8796093022207.0)", + "limits/memory", 0.0, (SCIP_Real)SCIP_MEM_NOLIMIT); + + AddSolverOption("lim:nodes", + "Maximal number of nodes to process (default: -1: no limit)", + "limits/nodes", -1, (int)SCIP_LONGINT_MAX); + + AddSolverOption("lim:restarts", + "Solving stops, if the given number of restarts was triggered (default: -1: no limit)", + "limits/restarts", -1, INT_MAX); + + AddSolverOption("lim:softtime softtime", + "Soft time limit which should be applied after first solution was found (default: -1.0: disabled)", + "limits/softtime", -1.0, SCIP_REAL_MAX); + + AddSolverOption("lim:solutions", + "Solving stops, if the given number of solutions were found (default: -1: no limit)", + "limits/solutions", -1, INT_MAX); + + AddSolverOption("lim:stallnodes", + "Solving stops, if the given number of nodes was processed since the last improvement of the primal solution value (default: -1: no limit)", + "limits/stallnodes", -1, (int)SCIP_LONGINT_MAX); + + AddSolverOption("lim:time timelim timelimit time_limit", + "Limit on solve time (in seconds; default: 1e+20).", + "limits/time", 0.0, 1e+20); + + AddSolverOption("lim:totalnodes", + "Maximal number of total nodes (incl. restarts) to process (default: -1: no limit)", + "limits/totalnodes", -1, (int)SCIP_LONGINT_MAX); + + ////////////////////////// LP ////////////////////////// + AddSolverOption("lp:pricing pricing", + "Pricing strategy:\n" + "\n.. value-table::", + "lp/pricing", values_pricing, "l"); + + AddSolverOption("lp:presolving", + "0/1: whether presolving of LP solver should be used" + "\n" + " | 0 - Presolving of LP solver should not be used\n" + " | 1 - Presolving of LP solver should be used (default).", + "lp/advanced/presolving", 0, 1); + + AddSolverOption("lp:threads", + "Number of threads used for solving the LP (default: 0: automatic)", + "lp/advanced/threads", 0, 64); + + AddSolverOption("lp:alwaysgetduals alwaysgetfarkasduals alwaysgetduals", + "0/1: whether the Farkas duals should always be collected when an LP is found to be infeasible" + "\n" + " | 0 - The Farkas duals should not always be collected when an LP is found to be infeasible (default)\n" + " | 1 - The Farkas duals should always be collected when an LP is found to be infeasible.", + "lp/alwaysgetduals", 0, 1); + + AddSolverOption("lp:solvedepth", + "Maximal depth for solving LP at the nodes (default: -1: no depth limit)", + "lp/solvedepth", -1, SCIP_MAXTREEDEPTH); + + AddSolverOption("lp:solvefreq", + "Frequency for solving LP at the nodes (-1: never; 0: only root LP; default: 1)", + "lp/solvefreq", -1, SCIP_MAXTREEDEPTH); + + ///////////////////////// MISC ///////////////////////// + AddSolverOption("misc:allowstrongdualreds allowstrongdualreds", + "0/1: whether strong dual reductions should be allowed in propagation and presolving" + "\n" + " | 0 - Strong dual reductions should not be allowed in propagation and presolving\n" + " | 1 - Strong dual reductions should be allowed in propagation and presolving (default).", + "misc/allowstrongdualreds", 0, 1); + + AddSolverOption("misc:allowweakdualreds allowweakdualreds", + "0/1: whether weak dual reductions should be allowed in propagation and presolving" + "\n" + " | 0 - Weak dual reductions should not be allowed in propagation and presolving\n" + " | 1 - Weak dual reductions should be allowed in propagation and presolving (default).", + "misc/allowweakdualreds", 0, 1); + + AddSolverOption("misc:scaleobj scaleobj", + "0/1: whether the objective function should be scaled so that it is always integer" + "\n" + " | 0 - The objective function should not be scaled so that it is always integer\n" + " | 1 - The objective function should be scaled so that it is always integer (default).", + "misc/scaleobj", 0, 1); + + ////////////////////////// NLP ///////////////////////// + AddSolverOption("nlp:disable", + "0/1: whether the NLP relaxation should be always disabled (also for NLPs/MINLPs)" + "\n" + " | 0 - NLP relaxation should not be always disabled (default)\n" + " | 1 - NLP relaxation should be always disabled.", + "nlp/disable", 0, 1); + + ////////////////////// NUMERICS //////////////////////// + AddSolverOption("num:checkfeastolfac checkfeastolfac", + "Feasibility tolerance factor; for checking the feasibility of the best solution (default: 1.0)", + "numerics/checkfeastolfac", 0.0, SCIP_REAL_MAX); + + AddSolverOption("num:dualfeastol dualfeastol", + "Feasibility tolerance for reduced costs in LP solution (default: 1e-07)", + "numerics/dualfeastol", SCIP_MINEPSILON*1e+03, SCIP_MAXEPSILON); + + AddSolverOption("num:epsilon epsilon", + "Absolute values smaller than this are considered zero (default: 1e-09)", + "numerics/epsilon", SCIP_MINEPSILON, SCIP_MAXEPSILON); + + AddSolverOption("num:feastol feastol", + "Feasibility tolerance for constraints (default: 1e-06)", + "numerics/feastol", SCIP_MINEPSILON*1e+03, SCIP_MAXEPSILON); + + AddSolverOption("num:infinity infinity", + "Values larger than this are considered infinity (default: 1e+20)", + "numerics/infinity", 1e+10, SCIP_INVALID/10.0); + + AddSolverOption("num:lpfeastolfactor lpfeastolfactor", + "Factor w.r.t. primal feasibility tolerance that determines default (and maximal) primal feasibility tolerance of LP solver (default: 1.0)", + "numerics/lpfeastolfactor", 1e-6, 1.0); + + AddSolverOption("num:sumepsilon sumepsilon", + "Absolute values of sums smaller than this are considered (default: 1e-06)", + "numerics/sumepsilon", SCIP_MINEPSILON*1e+03, SCIP_MAXEPSILON); + + //////////////////// NODESELECTION ///////////////////// + AddSolverOption("nod:childsel", + "Child selection rule:\n" + "\n.. value-table::\n", "nodeselection/childsel", childsel, "h"); + + ////////////////////// PARALLEL //////////////////////// + AddSolverOption("par:maxnthreads maxnthreads", + "Maximum number of threads used during parallel solve (default: 8)", + "parallel/maxnthreads", 0, 64); + + AddSolverOption("par:minnthreads minnthreads", + "Minimum number of threads used during parallel solve (default: 1)", + "parallel/minnthreads", 0, 64); + + AddSolverOption("par:mode mode", + "0/1: Parallel optimisation mode" + "\n" + " | 0 - Opportunistic\n" + " | 1 - Deterministic (default)", + "parallel/mode", 0, 1); + + ////////////////////// PRESOLVE //////////////////////// + AddSolverOption("pre:abortfac abortfac", + "Abort presolve, if at most this fraction of the problem was changed in last presolve round (default: 0.0008)", + "presolving/advanced/abortfac", 0.0, 1.0); + + AddSolverOption("pre:clqtablefac clqtablefac", + "Limit on number of entries in clique table relative to number of problem nonzeros (default: 2.0)", + "presolving/advanced/clqtablefac", 0.0, SCIP_REAL_MAX); + + AddSolverOption("pre:donotaggr donotaggr", + "0/1: whether aggregation of variables should be forbidden" + "\n" + " | 0 - Aggregation of variables should not be forbidden (default)\n" + " | 1 - Aggregation of variables should be forbidden.", + "presolving/advanced/donotaggr", 0, 1); + + AddSolverOption("pre:donotmultaggr donotmultaggr", + "0/1: whether multi-aggregation of variables should be forbidden" + "\n" + " | 0 - Multi-aggregation of variables should not be forbidden (default)\n" + " | 1 - Multi-aggregation of variables should be forbidden.", + "presolving/advanced/donotmultaggr", 0, 1); + + AddSolverOption("pre:immrestartfac immrestartfac", + "Fraction of integer variables that were fixed in the root node triggering an immediate restart with preprocessing (default: 0.1)", + "presolving/advanced/immrestartfac", 0.0, 1.0); + + AddSolverOption("pre:restartfac restartfac", + "Fraction of integer variables that were fixed in the root node triggering a restart with preprocessing after root node evaluation (default: 0.025)", + "presolving/advanced/restartfac", 0.0, 1.0); + + AddSolverOption("pre:restartminred restartminred", + "Minimal fraction of integer variables removed after restart to allow for an additional restart (default: 0.1)", + "presolving/advanced/restartminred", 0.0, 1.0); + + AddSolverOption("pre:subrestartfac subrestartfac", + "Fraction of integer variables that were globally fixed during the solving process triggering a restart with preprocessing (default: 1.0)", + "presolving/advanced/subrestartfac", 0.0, 1.0); + + AddSolverOption("pre:maxrestarts", + "Maximal number of restarts (default: -1: unlimited)", + "presolving/maxrestarts", -1, INT_MAX); + + AddSolverOption("pre:maxrounds", + "Maximal number of presolving rounds (default: -1: unlimited; 0: off)", + "presolving/maxrounds", -1, INT_MAX); + + AddStoredOption("pre:settings", + "0/1/2/3: sets presolvings settings" + "\n" + " | 0 - Sets presolvings default (default)\n" + " | 1 - Sets presolvings aggressive\n" + " | 2 - Sets presolvings fast\n" + " | 3 - Sets presolvings off.", + storedOptions_.heuristics_, 0, 3); + + ////////////////////// PROPAGATE /////////////////////// + AddSolverOption("pro:abortoncutoff", + "0/1: whether propagation should be aborted immediately (setting this to 0 could help conflict analysis to produce more conflict constraints)" + "\n" + " | 0 - Propagation should not be aborted immediately\n" + " | 1 - Propagation should be aborted immediately (default).", + "propagating/abortoncutoff", 0, 1); + + AddSolverOption("pro:maxrounds", + "Maximal number of propagation rounds per node (-1: unlimited; 0: off; default: 100)", + "propagating/maxrounds", -1, INT_MAX); + + AddSolverOption("pro:maxroundsroot", + "Maximal number of propagation rounds in the root node (-1: unlimited; 0: off; default: 1000)", + "propagating/maxroundsroot", -1, INT_MAX); + + //////////////////// RANDOMIZATION ///////////////////// + AddSolverOption("ran:permuteconss permuteconss", + "0/1: whether the order of constraints should be permuted (depends on permutationseed)? " + "\n" + " | 0 - Order of constraints should not be permuted\n" + " | 1 - Order of constraints should be permuted (default).", + "randomization/advanced/permuteconss", 0, 1); + + AddSolverOption("ran:permutevars permutevars", + "0/1: whether the order of variables should be permuted (depends on permutationseed)? " + "\n" + " | 0 - Order of variables should not be permuted (default)\n" + " | 1 - Order of variables should be permuted.", + "randomization/advanced/permutevars", 0, 1); + + AddSolverOption("ran:lpseed lpseed", + "Random seed for LP solver, e.g. for perturbations in the simplex (default: 0: LP default)", + "randomization/lpseed", 0, INT_MAX); + + AddSolverOption("ran:permutationseed permutationseed", + "Seed value for permuting the problem after reading/transformation (default: 0: no permutation) ", + "randomization/permutationseed", 0, INT_MAX); + + AddSolverOption("ran:randomseedshift randomseedshift", + "Global shift of all random seeds in the plugins and the LP random seed (default: 0) ", + "randomization/randomseedshift", 0, INT_MAX); + + ///////////////////////// TREE ///////////////////////// + AddSolverOption("est:method", + "Tree size estimation method:\n" + "\n.. value-table::\n", "estimation/method", estimation_method, "w"); + + AddSolverOption("est:completiontype", + "Approximation of search tree completion:\n" + "\n.. value-table::\n", "estimation/completiontype", estimation_completion, "a"); +} + + +double ScipBackend::MIPGap() { + return SCIPgetGap(getSCIP()) ScipBackend::VarStatii() { + + std::vector vars(NumVars()); + /* + SCIP_GetBasis(lp(), vars.data(), NULL); + for (auto& s : vars) { + switch (s) { + case SCIP_BASIS_BASIC: + s = (int)BasicStatus::bas; + break; + case SCIP_BASIS_LOWER: + s = (int)BasicStatus::low; + break; + case SCIP_BASIS_UPPER: + s = (int)BasicStatus::upp; + break; + case SCIP_BASIS_SUPERBASIC: + s = (int)BasicStatus::sup; + break; + case SCIP_BASIS_FIXED: + s = (int)BasicStatus::equ; + break; + default: + MP_RAISE(fmt::format("Unknown Scip VBasis value: {}", s)); + } + } + */ + return vars; +} + +ArrayRef ScipBackend::ConStatii() { + + std::vector cons(NumLinCons()); + /* + SCIP_GetBasis(lp(), NULL, cons.data()); + for (auto& s : cons) { + switch (s) { + case SCIP_BASIS_BASIC: + s = (int)BasicStatus::bas; + break; + case SCIP_BASIS_LOWER: + s = (int)BasicStatus::low; + break; + case SCIP_BASIS_UPPER: + s = (int)BasicStatus::upp; + break; + case SCIP_BASIS_SUPERBASIC: + s = (int)BasicStatus::sup; + break; + case SCIP_BASIS_FIXED: + s = (int)BasicStatus::equ; + break; + default: + MP_RAISE(fmt::format("Unknown Scip VBasis value: {}", s)); + } + }*/ + return cons; +} + +void ScipBackend::VarStatii(ArrayRef vst) { + int index[1]; + std::vector stt(vst.data(), vst.data() + vst.size()); + /* + for (auto j = stt.size(); j--; ) { + auto& s = stt[j]; + switch ((BasicStatus)s) { + case BasicStatus::bas: + s = SCIP_BASIS_BASIC; + break; + case BasicStatus::low: + s = SCIP_BASIS_LOWER; + break; + case BasicStatus::equ: + s = SCIP_BASIS_FIXED; + break; + case BasicStatus::upp: + s = SCIP_BASIS_UPPER; + break; + case BasicStatus::sup: + case BasicStatus::btw: + s = SCIP_BASIS_SUPERBASIC; + break; + case BasicStatus::none: + /// 'none' is assigned to new variables. Compute low/upp/sup: + /// Depending on where 0.0 is between bounds + double lb, ub; + index[0] = (int)j; + if(!SCIP_GetColInfo(lp(), SCIP_DBLINFO_LB, 1, index, &lb) && + !SCIP_GetColInfo(lp(), SCIP_DBLINFO_UB, 1, index, &ub)) + { + if (lb >= -1e-6) + s = -1; + else if (ub <= 1e-6) + s = -2; + else + s = -3; // or, leave at 0? + } + break; + default: + MP_RAISE(fmt::format("Unknown AMPL var status value: {}", s)); + } + } + SCIP_SetBasis(lp(), stt.data(), NULL); + */ +} + +void ScipBackend::ConStatii(ArrayRef cst) { + /* + std::vector stt(cst.data(), cst.data() + cst.size()); + for (auto& s : stt) { + switch ((BasicStatus)s) { + case BasicStatus::bas: + s = SCIP_BASIS_BASIC; + break; + case BasicStatus::none: // for 'none', which is the status + case BasicStatus::upp: // assigned to new rows, it seems good to guess + case BasicStatus::sup: // a valid status. + case BasicStatus::low: // + case BasicStatus::equ: // For active constraints, it is usually 'sup'. + case BasicStatus::btw: // We could compute slack to decide though. + s = SCIP_BASIS_SUPERBASIC; + break; + default: + MP_RAISE(fmt::format("Unknown AMPL con status value: {}", s)); + } + } + SCIP_SetBasis(lp(), NULL, stt.data()); + */ +} + +SolutionBasis ScipBackend::GetBasis() { + std::vector varstt = VarStatii(); + std::vector constt = ConStatii(); + if (varstt.size() && constt.size()) { + auto mv = GetValuePresolver().PostsolveBasis( + { std::move(varstt), + {{{ CG_Linear, std::move(constt) }}} }); + varstt = mv.GetVarValues()(); + constt = mv.GetConValues()(); + assert(varstt.size()); + } + return { std::move(varstt), std::move(constt) }; +} + +void ScipBackend::SetBasis(SolutionBasis basis) { + auto mv = GetValuePresolver().PresolveBasis( + { basis.varstt, basis.constt }); + auto varstt = mv.GetVarValues()(); + auto constt = mv.GetConValues()(CG_Linear); + assert(varstt.size()); + assert(constt.size()); + VarStatii(varstt); + ConStatii(constt); +} + +void ScipBackend::AddPrimalDualStart(Solution sol) +{ + auto mv = GetValuePresolver().PresolveSolution( + { sol.primal, sol.dual } ); + auto x0 = mv.GetVarValues()(); + auto pi0 = mv.GetConValues()(CG_Linear); + SCIP_SOL* solution; + SCIP_Bool keep; + SCIP_CCALL( SCIPcreateSol(getSCIP(), &solution, NULL) ); + + SCIP_CCALL( SCIPsetSolVals(getSCIP(), solution, getPROBDATA()->nvars, getPROBDATA()->vars, x0.data()) ); + + SCIP_CCALL( SCIPaddSolFree(getSCIP(), &solution, &keep) ); +} + +void ScipBackend::AddMIPStart(ArrayRef x0, ArrayRef sparsity) { + SCIP_SOL* solution; + SCIP_Bool keep; + SCIP_CCALL( SCIPcreateSol(getSCIP(), &solution, NULL) ); + + SCIP_CCALL( SCIPsetSolVals(getSCIP(), solution, getPROBDATA()->nvars, getPROBDATA()->vars, (double*)x0.data()) ); + + SCIP_CCALL( SCIPaddSolFree(getSCIP(), &solution, &keep) ); +} + + +} // namespace mp + + +// AMPLs +void* AMPLSOpenScip( + const char* slv_opt, CCallbacks cb = {}) { + return AMPLS__internal__Open(std::unique_ptr{new mp::ScipBackend()}, + slv_opt, cb); +} + +void AMPLSCloseScip(AMPLS_MP_Solver* slv) { + AMPLS__internal__Close(slv); +} + +void* GetScipmodel(AMPLS_MP_Solver* slv) { + return + dynamic_cast(AMPLSGetBackend(slv))->getSCIP(); +} diff --git a/solvers/scip/scipbackend.h b/solvers/scip/scipbackend.h new file mode 100644 index 000000000..bedcf3599 --- /dev/null +++ b/solvers/scip/scipbackend.h @@ -0,0 +1,168 @@ +#ifndef MP_SCIP_BACKEND_H_ +#define MP_SCIP_BACKEND_H_ + +#include + +#include "mp/backend-mip.h" +#include "mp/flat/backend_flat.h" +#include "scipcommon.h" + +namespace mp { + +class ScipBackend : + public FlatBackend< MIPBackend >, + public ScipCommon +{ + using BaseBackend = FlatBackend< MIPBackend >; + + //////////////////// [[ The public interface ]] ////////////////////// +public: + ScipBackend(); + ~ScipBackend(); + + /// Prefix used for the _options environment variable + static const char* GetAMPLSolverName() { return "scip"; } + + /// AMPL driver name displayed in messages + static const char* GetAMPLSolverLongName() { return "AMPL-SCIP"; } + /// Solver name displayed in messages + static const char* GetSolverName() { return "SCIP"; } + /// Version displayed with -v + std::string GetSolverVersion(); + + /// Name for diagnostic messages + static const char* GetBackendName(); + static const char* GetBackendLongName() { return nullptr; } + + /// Init custom driver options, such as outlev, writeprob + void InitCustomOptions() override; + /// Chance for the Backend to init solver environment, etc. + void InitOptionParsing() override { } + /// Chance to consider options immediately (open cloud, etc) + void FinishOptionParsing() override; + + + + //////////////////////////////////////////////////////////// + /////////////// OPTIONAL STANDARD FEATURES ///////////////// + //////////////////////////////////////////////////////////// + // Use this section to declare and implement some standard features + // that may or may not need additional functions. + USING_STD_FEATURES; + + /** + * MULTISOL support + * No API, see ReportIntermediateSolution() +**/ + ALLOW_STD_FEATURE(MULTISOL, false) + + /** + * Get/Set AMPL var/con statii + **/ + ALLOW_STD_FEATURE(BASIS, false) + // TODO If getting/setting a basis is supported, implement the + // accessor and the setter below + SolutionBasis GetBasis() override; + void SetBasis(SolutionBasis) override; + /** + * General warm start, e.g., + * set primal/dual initial guesses for continuous case + **/ + ALLOW_STD_FEATURE( WARMSTART, false ) + void AddPrimalDualStart(Solution sol) override; + /** + * MIP warm start + **/ + // If MIP warm start is supported, implement the function below + // to set a non-presolved starting solution + ALLOW_STD_FEATURE(MIPSTART, false) + void AddMIPStart(ArrayRef x0, + ArrayRef sparsity) override; + + + /** + * Get MIP Gap + **/ + // return MIP gap + // (adds option mip:return_gap) + ALLOW_STD_FEATURE(RETURN_MIP_GAP, true) + double MIPGap() override; + double MIPGapAbs() override; + /** + * Get MIP dual bound + **/ + // return the best dual bound value + // (adds option mip:bestbound) + ALLOW_STD_FEATURE(RETURN_BEST_DUAL_BOUND, true) + double BestDualBound() override; + + /////////////////////////// Model attributes ///////////////////////// + bool IsMIP() const override; + bool IsQCP() const override; + + //////////////////////////// SOLVING /////////////////////////////// + + /// Note the interrupt notifier + void SetInterrupter(mp::Interrupter* inter) override; + +public: // public for static polymorphism + /// Solve, no model modification any more (such as feasrelax). + /// Can report intermediate results via HandleFeasibleSolution() during this, + /// otherwise/finally via ReportResults() + void Solve() override; + + /// Default impl of GetObjValues() + ArrayRef GetObjectiveValues() override + { return std::vector{ObjectiveValue()}; } + + + //////////////////// [[ Implementation details ]] ////////////////////// + /////////////////////////////////////////////////////////////////////////////// +protected: + void ExportModel(const std::string& file); + + double ObjectiveValue() const; + + /// Solution values. The vectors are emptied if not available + ArrayRef PrimalSolution() override; + pre::ValueMapDbl DualSolution() override; + ArrayRef DualSolution_LP(); + + void WindupSCIPSolve(); + + void ReportResults() override; + void ReportSCIPResults(); + + void ReportSCIPPool(); + + std::vector getPoolSolution(int i); + double getPoolObjective(int i); + + /// Solution attributes + double NodeCount() const; + double SimplexIterations() const; + int BarrierIterations() const; + + std::pair ConvertSCIPStatus(); + void AddSCIPMessages(); + + ArrayRef VarStatii(); + ArrayRef ConStatii(); + void VarStatii(ArrayRef); + void ConStatii(ArrayRef); + + +private: + /// These options are stored in the class + struct Options { + std::string exportFile_, logFile_, paramRead_; + int concurrent_, heuristics_, cuts_, presolvings_ = 0; + }; + Options storedOptions_; + + +}; + +} // namespace mp + +#endif // MP_SCIP_BACKEND_H_ diff --git a/solvers/scip/scipcommon.cc b/solvers/scip/scipcommon.cc new file mode 100644 index 000000000..fa7a97f18 --- /dev/null +++ b/solvers/scip/scipcommon.cc @@ -0,0 +1,152 @@ +#include "mp/format.h" +#include "scipcommon.h" + +static +SCIP_DECL_PROBDELORIG(probdataDelOrigNl) +{ + int i; + + assert((*probdata)->vars != NULL || (*probdata)->nvars == 0); + assert((*probdata)->linconss != NULL || (*probdata)->nlinconss == 0); + + for( i = 0; i < (*probdata)->nlinconss; ++i ) + { + SCIP_CALL( SCIPreleaseCons(scip, &(*probdata)->linconss[i]) ); + } + SCIPfreeBlockMemoryArray(scip, &(*probdata)->linconss, (*probdata)->nlinconss); + + for( i = 0; i < (*probdata)->nvars; ++i ) + { + SCIP_CCALL( SCIPreleaseVar(scip, &(*probdata)->vars[i]) ); + } + SCIPfreeBlockMemoryArray(scip, &(*probdata)->vars, (*probdata)->nvars); + + SCIPfreeMemory(scip, probdata); + + return SCIP_OKAY; +} + +namespace mp { + +void ScipCommon::OpenSolver() { + int status = 0; + SCIP* scip = NULL; + SCIP_PROBDATA* probdata = NULL; + + // initialize SCIP + status = SCIPcreate(&scip); + setSCIP(scip); // Assign it + + // include default SCIP plugins + SCIP_CCALL( SCIPincludeDefaultPlugins(scip) ); + + // initialize empty SCIP problem + SCIP_CCALL( SCIPallocClearMemory(scip, &probdata) ); + SCIP_CCALL( SCIPcreateProb(scip, "", probdataDelOrigNl, NULL, NULL, NULL, NULL, NULL, probdata) ); + setPROBDATA(probdata); + + if (status != 1) + throw std::runtime_error( fmt::format( + "Failed to create problem, error code {}.", status ) ); + SetSolverOption("display/verblevel", 0); +} + +void ScipCommon::CloseSolver() { + SCIP* scip = getSCIP(); + + // free SCIP + SCIP_CCALL( SCIPfree(&scip) ); +} + +int ScipCommon::NumLinCons() const { + return getPROBDATA()->nlinconss; +} + +int ScipCommon::NumVars() const { + return getPROBDATA()->nvars; +} + +int ScipCommon::NumObjs() const { + return 1; +} + +int ScipCommon::NumQPCons() const { + // TODO Get number of quadratic constraints using solver API + // return getIntAttr(SCIP_INTATTR_QCONSTRS); + return 0; +} + +int ScipCommon::NumSOSCons() const { + // TODO Get number of SOS constraints using solver API + // return getIntAttr(SCIP_INTATTR_SOSS); + return 0; +} + +int ScipCommon::NumIndicatorCons() const { + // TODO Get number of indicator constraints using solver API + // return getIntAttr(SCIP_INTATTR_INDICATORS); + return 0; +} + +void ScipCommon::GetSolverOption(const char* key, int &value) const { + if (SCIPparamGetType(SCIPgetParam(getSCIP(), key))==SCIP_PARAMTYPE_BOOL) { + SCIP_Bool buffer; + SCIP_CCALL( SCIPgetBoolParam(getSCIP(), key, &buffer) ); + value = (int)buffer; + } + else + SCIP_CCALL( SCIPgetIntParam(getSCIP(), key, &value) ); +} + +void ScipCommon::SetSolverOption(const char* key, int value) { + if (SCIPparamGetType(SCIPgetParam(getSCIP(), key))==SCIP_PARAMTYPE_BOOL) + SCIP_CCALL( SCIPsetBoolParam(getSCIP(), key, value) ); + else + SCIP_CCALL( SCIPsetIntParam(getSCIP(), key, value) ); +} + +void ScipCommon::GetSolverOption(const char* key, double &value) const { + SCIP_CCALL( SCIPgetRealParam(getSCIP(), key, &value) ); +} + +void ScipCommon::SetSolverOption(const char* key, double value) { + SCIP_CCALL( SCIPsetRealParam(getSCIP(), key, value) ); +} + +void ScipCommon::GetSolverOption(const char* key, std::string &value) const { + char* buffer; + SCIP_CCALL( SCIPallocBufferArray(getSCIP(), &buffer, SCIP_MAXSTRLEN) ); + if (SCIPparamGetType(SCIPgetParam(getSCIP(), key))==SCIP_PARAMTYPE_CHAR) + SCIP_CCALL( SCIPgetCharParam(getSCIP(), key, buffer) ); + else + SCIP_CCALL( SCIPgetStringParam(getSCIP(), key, &buffer) ); + value = buffer; + SCIPfreeBufferArray(getSCIP(), &buffer); +} + +void ScipCommon::SetSolverOption(const char* key, const std::string& value) { + if (SCIPparamGetType(SCIPgetParam(getSCIP(), key))==SCIP_PARAMTYPE_CHAR) + SCIP_CCALL( SCIPsetCharParam(getSCIP(), key, value.c_str()[0]) ); + else + SCIP_CCALL( SCIPsetStringParam(getSCIP(), key, value.c_str()) ); +} + + +double ScipCommon::Infinity() const { + double inf; + GetSolverOption("numerics/infinity", inf); + return inf; +} + +double ScipCommon::MinusInfinity() { return -Infinity(); } + +bool ScipCommon::IsContinuous() { + // True iff scip has only continuous variables + for (int i = 0; i < getPROBDATA()->nvars; i++) { + if (SCIPvarGetType(getPROBDATA()->vars[i]) != SCIP_VARTYPE_CONTINUOUS) + return false; + } + return true; +} + +} // namespace mp diff --git a/solvers/scip/scipcommon.h b/solvers/scip/scipcommon.h new file mode 100644 index 000000000..0a0fbbe17 --- /dev/null +++ b/solvers/scip/scipcommon.h @@ -0,0 +1,89 @@ +#ifndef SCIPCOMMON_H +#define SCIPCOMMON_H + +#include + +#include "mp/backend-to-model-api.h" + +#include "scip/scip.h" +#include "scip/scipdefplugins.h" + +#include "mp/format.h" + +/// problem data stored in SCIP +struct SCIP_ProbData +{ + SCIP_VAR** vars; /**< variables in the order given by AMPL */ + int nvars; /**< number of variables */ + + SCIP_CONS** linconss; /**< linear constraints in the order given by AMPL */ + int i; /**< shows free slot of linear constraints */ + int nlinconss; /**< number of linear constraints */ +}; + +namespace mp { + +/// Information shared by both +/// `ScipBackend` and `ScipModelAPI` +struct ScipCommonInfo { + SCIP* getSCIP() const { return scip_; } + void setSCIP(SCIP* scip) { scip_ = scip; } + + SCIP_PROBDATA* getPROBDATA() const { return probdata_; } + void setPROBDATA(SCIP_PROBDATA* probdata) { probdata_ = probdata; } + +private: + SCIP* scip_ = NULL; + SCIP_PROBDATA* probdata_; +}; + + +/// Common API for Scip classes +class ScipCommon : + public Backend2ModelAPIConnector { +public: + /// These methods access Scip options. Used by AddSolverOption() + void GetSolverOption(const char* key, int& value) const; + void SetSolverOption(const char* key, int value); + void GetSolverOption(const char* key, double& value) const; + void SetSolverOption(const char* key, double value); + void GetSolverOption(const char* key, std::string& value) const; + void SetSolverOption(const char* key, const std::string& value); + + /// SCIP own infinity + double Infinity() const; + double MinusInfinity(); + + bool IsContinuous(); + +protected: + void OpenSolver(); + void CloseSolver(); + + int NumLinCons() const; + int NumVars() const; + int NumObjs() const; + int NumQPCons() const; + int NumSOSCons() const; + int NumIndicatorCons() const; + +protected: + // TODO if desirable, provide function to create the solver's environment + // with own license + // int (*createEnv) (solver_env**) = nullptr; + +}; + + +/// Convenience macro +// TODO This macro is useful to automatically throw an error if a function in the +// solver API does not return a valid errorcode. In this mock driver, we define it +// ourselves, normally this constant would be defined in the solver's API. +#define SCIP_RETCODE_OK 1 +#define SCIP_CCALL( call ) do { if (int e = (call) != SCIP_RETCODE_OK) \ + throw std::runtime_error( \ + fmt::format(" Call failed: '{}' with code {}", #call, e )); } while (0) + +} // namespace mp + +#endif // SCIPCOMMON_H diff --git a/solvers/scip/scipmodelapi.cc b/solvers/scip/scipmodelapi.cc new file mode 100644 index 000000000..2cf118efb --- /dev/null +++ b/solvers/scip/scipmodelapi.cc @@ -0,0 +1,455 @@ +#include "scipmodelapi.h" + +namespace mp { + +void ScipModelAPI::InitProblemModificationPhase(const FlatModelInfo* flat_model_info) { + // Allocate storage if needed: + int n_linear_cons = flat_model_info->GetNumberOfConstraintsOfGroup(CG_Linear); + getPROBDATA()->nlinconss = n_linear_cons; + getPROBDATA()->i = 0; + SCIP_CCALL( SCIPallocBlockMemoryArray(getSCIP(), &getPROBDATA()->linconss, getPROBDATA()->nlinconss) ); +} + +void ScipModelAPI::AddVariables(const VarArrayDef& v) { + assert( SCIPgetStage(getSCIP()) == SCIP_STAGE_PROBLEM ); + + getPROBDATA()->nvars = v.size(); + SCIP_CCALL( SCIPallocBlockMemoryArray(getSCIP(), &getPROBDATA()->vars, getPROBDATA()->nvars) ); + + for (int i = 0; i < v.size(); i++) { + SCIP_VAR* var = NULL; + double lb = std::isinf(v.plb()[i]) ? MinusInfinity() : v.plb()[i]; + double ub = std::isinf(v.pub()[i]) ? Infinity() : v.pub()[i]; + double objcoef = 0.0; + SCIP_VARTYPE vartype; + if (v.ptype()[i] == var::Type::INTEGER && lb == 0.0 && ub == 1.0) + vartype = SCIP_VARTYPE_BINARY; + else if (v.ptype()[i] == var::Type::INTEGER) + vartype = SCIP_VARTYPE_INTEGER; + else + vartype = SCIP_VARTYPE_CONTINUOUS; + //const char* name = v.pnames()[i]; + SCIP_CCALL( SCIPcreateVarBasic(getSCIP(), &var, NULL, lb, ub, objcoef, vartype) ); + SCIP_CCALL( SCIPaddVar(getSCIP(), var) ); + getPROBDATA()->vars[i] = var; + } +} + +void ScipModelAPI::SetLinearObjective( int iobj, const LinearObjective& lo ) { + if (iobj<1) { + SCIP_CCALL( SCIPsetObjsense(getSCIP(), + obj::Type::MAX==lo.obj_sense() ? SCIP_OBJSENSE_MAXIMIZE : SCIP_OBJSENSE_MINIMIZE) ); + SCIP_VAR** vars = getPROBDATA()->vars; + for (int i = 0; i < lo.num_terms(); i++) { + SCIP_CCALL( SCIPchgVarObj(getSCIP(), vars[lo.vars()[i]], lo.coefs()[i]) ); + } + } + else { + throw std::runtime_error("Multiple linear objectives not supported"); + } +} + + +void ScipModelAPI::SetQuadraticObjective(int iobj, const QuadraticObjective& qo) { + throw std::runtime_error("Quadratic objective not supported"); +} + + +void ScipModelAPI::linearHelper(const int* pvars, const double* pcoefs, const size_t size, const char* name, const double lb, const double ub) { + SCIP_VAR** vars = NULL; + SCIP_CCALL( SCIPallocBufferArray(getSCIP(), &vars, size) ); + for (size_t i = 0; i < size; i++) + vars[i] = getPROBDATA()->vars[pvars[i]]; + + SCIP_CONS* cons; + SCIP_CCALL( SCIPcreateConsBasicLinear(getSCIP(), &cons, name, (int)size, vars, (double*)pcoefs, lb, ub) ); + SCIP_CCALL( SCIPaddCons(getSCIP(), cons) ); + getPROBDATA()->linconss[getPROBDATA()->i] = cons; + getPROBDATA()->i++; + + SCIPfreeBufferArray(getSCIP(), &vars); +} + +void ScipModelAPI::AddConstraint(const LinConRange& lc) { + linearHelper(lc.pvars(), lc.pcoefs(), lc.size(), lc.GetName(), lc.lb(), lc.ub()); +} + +void ScipModelAPI::AddConstraint(const LinConLE& lc) { + linearHelper(lc.pvars(), lc.pcoefs(), lc.size(), lc.GetName(), MinusInfinity(), lc.ub()); +} + +void ScipModelAPI::AddConstraint(const LinConEQ& lc) { + linearHelper(lc.pvars(), lc.pcoefs(), lc.size(), lc.GetName(), lc.ub(), lc.ub()); +} + +void ScipModelAPI::AddConstraint(const LinConGE& lc) { + linearHelper(lc.pvars(), lc.pcoefs(), lc.size(), lc.GetName(), lc.lb(), Infinity()); +} + + +void ScipModelAPI::AddConstraint(const AbsConstraint &absc) { + SCIP_VAR* x = getPROBDATA()->vars[absc.GetArguments()[0]]; + SCIP_VAR* res = getPROBDATA()->vars[absc.GetResultVar()]; + + assert(x != NULL); + assert(res != NULL); + + SCIP_EXPR* terms[2]; + SCIP_EXPR* xexpr; + SCIP_EXPR* sumexpr; + SCIP_Real coefs[2] = {1.0, -1.0}; + + SCIP_CCALL( SCIPcreateExprVar(getSCIP(), &xexpr, x, NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprVar(getSCIP(), &terms[0], res, NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprAbs(getSCIP(), &terms[1], xexpr, NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprSum(getSCIP(), &sumexpr, 2, terms, coefs, 0.0, NULL, NULL) ); + + SCIP_CONS* cons; + SCIP_CCALL( SCIPcreateConsBasicNonlinear(getSCIP(), &cons, absc.GetName(), sumexpr, 0.0, 0.0) ); + SCIP_CCALL( SCIPaddCons(getSCIP(), cons) ); + SCIP_CCALL( SCIPreleaseCons(getSCIP(), &cons) ); + + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &sumexpr) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &terms[1]) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &terms[0]) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &xexpr) ); +} + +void ScipModelAPI::AddConstraint(const AndConstraint &cc) { + SCIP_VAR** vars = NULL; + SCIP_CCALL( SCIPallocBufferArray(getSCIP(), &vars, cc.GetArguments().size()) ); + for (size_t i = 0; i < cc.GetArguments().size(); i++) { + vars[i] = getPROBDATA()->vars[cc.GetArguments()[i]]; + } + + SCIP_CONS* cons; + SCIP_CCALL( SCIPcreateConsBasicAnd(getSCIP(), &cons, cc.GetName(), getPROBDATA()->vars[cc.GetResultVar()], + (int)cc.GetArguments().size(), vars) ); + SCIP_CCALL( SCIPaddCons(getSCIP(), cons) ); + SCIP_CCALL( SCIPreleaseCons(getSCIP(), &cons) ); + + SCIPfreeBufferArray(getSCIP(), &vars); +} + +void ScipModelAPI::AddConstraint(const OrConstraint &dc) { + SCIP_VAR** vars = NULL; + SCIP_CCALL( SCIPallocBufferArray(getSCIP(), &vars, dc.GetArguments().size()) ); + for (size_t i = 0; i < dc.GetArguments().size(); i++) { + vars[i] = getPROBDATA()->vars[dc.GetArguments()[i]]; + } + + SCIP_CONS* cons; + SCIP_CCALL( SCIPcreateConsBasicOr(getSCIP(), &cons, dc.GetName(), getPROBDATA()->vars[dc.GetResultVar()], + (int)dc.GetArguments().size(), vars) ); + SCIP_CCALL( SCIPaddCons(getSCIP(), cons) ); + SCIP_CCALL( SCIPreleaseCons(getSCIP(), &cons) ); + + SCIPfreeBufferArray(getSCIP(), &vars); +} + + +void ScipModelAPI::helpIndicatorLin(const int* pvars, const double* pcoefs, const size_t size, const char* name, const double rhs, const int binary_var, const int binary_value, SCIP_Bool lessthanineq) { + SCIP_VAR** vars = NULL; + SCIP_CCALL( SCIPallocBufferArray(getSCIP(), &vars, size) ); + for (size_t i = 0; i < size; i++) + vars[i] = getPROBDATA()->vars[pvars[i]]; + SCIP_CONS* cons; + SCIP_CCALL( SCIPcreateConsIndicatorGeneric(getSCIP(), &cons, name, getPROBDATA()->vars[binary_var], + (int)size, vars, (double*)pcoefs, rhs, binary_value, lessthanineq, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE) ); + SCIP_CCALL( SCIPaddCons(getSCIP(), cons) ); + SCIP_CCALL( SCIPreleaseCons(getSCIP(), &cons) ); + + SCIPfreeBufferArray(getSCIP(), &vars); +} + +void ScipModelAPI::AddConstraint(const IndicatorConstraintLinLE &ic) { + helpIndicatorLin(ic.get_constraint().pvars(), ic.get_constraint().pcoefs(), ic.get_constraint().size(), ic.GetName(), ic.get_constraint().rhs(), ic.get_binary_var(), ic.get_binary_value(), TRUE); +} + +void ScipModelAPI::AddConstraint(const IndicatorConstraintLinEQ &ic) { + SCIP_VAR** vars = NULL; + SCIP_CCALL( SCIPallocBufferArray(getSCIP(), &vars, ic.get_constraint().size()) ); + for (size_t i = 0; i < ic.get_constraint().size(); i++) { + vars[i] = getPROBDATA()->vars[ic.get_constraint().pvars()[i]]; + } + SCIP_CONS* cons1; + SCIP_CONS* cons2; + SCIP_CCALL( SCIPcreateConsIndicatorGeneric(getSCIP(), &cons1, ic.GetName(),getPROBDATA()->vars[ic.get_binary_var()], + (int)ic.get_constraint().size(), vars, (double*)ic.get_constraint().pcoefs(), ic.get_constraint().rhs(), + ic.get_binary_value(), TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE) ); + SCIP_CCALL( SCIPcreateConsIndicatorGeneric(getSCIP(), &cons2, ic.GetName(),getPROBDATA()->vars[ic.get_binary_var()], + (int)ic.get_constraint().size(), vars, (double*)ic.get_constraint().pcoefs(), ic.get_constraint().rhs(), + ic.get_binary_value(), FALSE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE) ); + SCIP_CCALL( SCIPaddCons(getSCIP(), cons1) ); + SCIP_CCALL( SCIPaddCons(getSCIP(), cons2) ); + SCIP_CCALL( SCIPreleaseCons(getSCIP(), &cons1) ); + SCIP_CCALL( SCIPreleaseCons(getSCIP(), &cons2) ); + + SCIPfreeBufferArray(getSCIP(), &vars); +} + +void ScipModelAPI::AddConstraint(const IndicatorConstraintLinGE &ic) { + helpIndicatorLin(ic.get_constraint().pvars(), ic.get_constraint().pcoefs(), ic.get_constraint().size(), ic.GetName(), ic.get_constraint().rhs(), ic.get_binary_var(), ic.get_binary_value(), FALSE); +} + + +void ScipModelAPI::helpQuad(const char* name, const int* lt_pvars, const double* lt_pcoefs, const size_t lt_size, const int* qt_pvars1, const int* qt_pvars2, const double* qt_pcoefs, const size_t qt_size, const double lb, const double ub) { + // allocation of linterm vars + SCIP_VAR** linvars = NULL; + SCIP_CCALL( SCIPallocBufferArray(getSCIP(), &linvars, lt_size) ); + for (size_t i = 0; i < lt_size; i++) + linvars[i] = getPROBDATA()->vars[lt_pvars[i]]; + + // allocation of qterm1 vars + SCIP_VAR** qvars1 = NULL; + SCIP_CCALL( SCIPallocBufferArray(getSCIP(), &qvars1, qt_size) ); + for (size_t i = 0; i < qt_size; i++) + qvars1[i] = getPROBDATA()->vars[qt_pvars1[i]]; + + // allocation of qterm2 vars + SCIP_VAR** qvars2 = NULL; + SCIP_CCALL( SCIPallocBufferArray(getSCIP(), &qvars2, qt_size) ); + for (size_t i = 0; i < qt_size; i++) + qvars2[i] = getPROBDATA()->vars[qt_pvars2[i]]; + + + SCIP_CONS* cons; + SCIP_CCALL( SCIPcreateConsBasicQuadraticNonlinear(getSCIP(), &cons, name, (int)lt_size, linvars, (double*)lt_pcoefs, (int)qt_size, qvars1, qvars2, (double*)qt_pcoefs, lb, ub) ); + SCIP_CCALL( SCIPaddCons(getSCIP(), cons) ); + SCIP_CCALL( SCIPreleaseCons(getSCIP(), &cons) ); + + SCIPfreeBufferArray(getSCIP(), &linvars); + SCIPfreeBufferArray(getSCIP(), &qvars1); + SCIPfreeBufferArray(getSCIP(), &qvars2); +} + +void ScipModelAPI::AddConstraint(const QuadConRange& qc) { + const auto& lt = qc.GetLinTerms(); + const auto& qt = qc.GetQPTerms(); + + helpQuad(qc.GetName(), lt.pvars(), lt.pcoefs(), lt.size(), qt.pvars1(), qt.pvars2(), qt.pcoefs(), qt.size(), qc.lb(), qc.ub()); +} + +void ScipModelAPI::AddConstraint( const QuadConLE& qc ) { + const auto& lt = qc.GetLinTerms(); + const auto& qt = qc.GetQPTerms(); + + helpQuad(qc.GetName(), lt.pvars(), lt.pcoefs(), lt.size(), qt.pvars1(), qt.pvars2(), qt.pcoefs(), qt.size(), MinusInfinity(), qc.rhs()); +} + +void ScipModelAPI::AddConstraint( const QuadConEQ& qc ) { + const auto& lt = qc.GetLinTerms(); + const auto& qt = qc.GetQPTerms(); + + helpQuad(qc.GetName(), lt.pvars(), lt.pcoefs(), lt.size(), qt.pvars1(), qt.pvars2(), qt.pcoefs(), qt.size(), qc.rhs(), qc.rhs()); +} + +void ScipModelAPI::AddConstraint( const QuadConGE& qc ) { + const auto& lt = qc.GetLinTerms(); + const auto& qt = qc.GetQPTerms(); + + helpQuad(qc.GetName(), lt.pvars(), lt.pcoefs(), lt.size(), qt.pvars1(), qt.pvars2(), qt.pcoefs(), qt.size(), qc.lb(), Infinity()); +} + +/* +void ScipModelAPI::AddConstraint( const QuadraticConeConstraint& qc ) { + const auto& arg = qc.GetArguments(); + SCIP_VAR** vars = NULL; + SCIP_CCALL( SCIPallocBufferArray(getSCIP(), &vars, arg.size()) ); + for (size_t i = 0; i < arg.size(); i++) { + vars[i] = getPROBDATA()->vars[arg[i]]; + } + + SCIP_CONS* cons; + SCIP_CCALL( SCIPcreateConsBasicSOCNonlinear(getSCIP(), &cons, qc.GetName(), (int)arg.size(), vars, (SCIP_Real*)qc.GetParameters().data(), NULL, 0.0, getPROBDATA()->vars[qc.GetResultVar()], 1.0, 0.0) ); + SCIP_CCALL( SCIPaddCons(getSCIP(), cons) ); + SCIP_CCALL( SCIPreleaseCons(getSCIP(), &cons) ); + + SCIPfreeBufferArray(getSCIP(), &vars); +}*/ + + +void ScipModelAPI::AddConstraint(const SOS1Constraint& sos) { + SCIP_VAR** vars = NULL; + double* weights = NULL; + SCIP_CCALL( SCIPallocBufferArray(getSCIP(), &vars, sos.size()) ); + SCIP_CCALL( SCIPallocBufferArray(getSCIP(), &weights, sos.size()) ); + for (int i = 0; i < sos.size(); i++) { + vars[i] = getPROBDATA()->vars[sos.get_vars().data()[i]]; + weights[i] = sos.get_weights().data()[i]; + } + + SCIP_CONS* cons = NULL; + SCIP_CCALL( SCIPcreateConsBasicSOS1(getSCIP(), &cons, sos.GetName(), sos.size(), vars, weights) ); + SCIP_CCALL( SCIPaddCons(getSCIP(), cons) ); + SCIP_CCALL( SCIPreleaseCons(getSCIP(), &cons) ); + + SCIPfreeBufferArray(getSCIP(), &vars); + SCIPfreeBufferArray(getSCIP(), &weights); +} + +void ScipModelAPI::AddConstraint(const SOS2Constraint& sos) { + SCIP_VAR** vars = NULL; + double* weights = NULL; + SCIP_CCALL( SCIPallocBufferArray(getSCIP(), &vars, sos.size()) ); + SCIP_CCALL( SCIPallocBufferArray(getSCIP(), &weights, sos.size()) ); + for (int i = 0; i < sos.size(); i++) { + vars[i] = getPROBDATA()->vars[sos.get_vars().data()[i]]; + weights[i] = sos.get_weights().data()[i]; + } + + SCIP_CONS* cons = NULL; + SCIP_CCALL( SCIPcreateConsBasicSOS2(getSCIP(), &cons, sos.GetName(), sos.size(), vars, weights) ); + SCIP_CCALL( SCIPaddCons(getSCIP(), cons) ); + SCIP_CCALL( SCIPreleaseCons(getSCIP(), &cons) ); + + SCIPfreeBufferArray(getSCIP(), &vars); + SCIPfreeBufferArray(getSCIP(), &weights); +} + + +void ScipModelAPI::AddConstraint(const ExpConstraint &cc) { + SCIP_VAR* x = getPROBDATA()->vars[cc.GetArguments()[0]]; + SCIP_VAR* res = getPROBDATA()->vars[cc.GetResultVar()]; + + assert(x != NULL); + assert(res != NULL); + + SCIP_EXPR* terms[2]; + SCIP_EXPR* xexpr; + SCIP_EXPR* sumexpr; + SCIP_Real coefs[2] = {1.0, -1.0}; + + SCIP_CCALL( SCIPcreateExprVar(getSCIP(), &xexpr, x, NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprVar(getSCIP(), &terms[0], res, NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprExp(getSCIP(), &terms[1], xexpr, NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprSum(getSCIP(), &sumexpr, 2, terms, coefs, 0.0, NULL, NULL) ); + + SCIP_CONS* cons; + SCIP_CCALL( SCIPcreateConsBasicNonlinear(getSCIP(), &cons, cc.GetName(), sumexpr, 0.0, 0.0) ); + SCIP_CCALL( SCIPaddCons(getSCIP(), cons) ); + SCIP_CCALL( SCIPreleaseCons(getSCIP(), &cons) ); + + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &sumexpr) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &terms[1]) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &terms[0]) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &xexpr) ); +} + +void ScipModelAPI::AddConstraint(const LogConstraint &cc) { + SCIP_VAR* x = getPROBDATA()->vars[cc.GetArguments()[0]]; + SCIP_VAR* res = getPROBDATA()->vars[cc.GetResultVar()]; + + assert(x != NULL); + assert(res != NULL); + + SCIP_EXPR* terms[2]; + SCIP_EXPR* xexpr; + SCIP_EXPR* sumexpr; + SCIP_Real coefs[2] = {1.0, -1.0}; + + SCIP_CCALL( SCIPcreateExprVar(getSCIP(), &xexpr, x, NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprVar(getSCIP(), &terms[0], res, NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprLog(getSCIP(), &terms[1], xexpr, NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprSum(getSCIP(), &sumexpr, 2, terms, coefs, 0.0, NULL, NULL) ); + + SCIP_CONS* cons; + SCIP_CCALL( SCIPcreateConsBasicNonlinear(getSCIP(), &cons, cc.GetName(), sumexpr, 0.0, 0.0) ); + SCIP_CCALL( SCIPaddCons(getSCIP(), cons) ); + SCIP_CCALL( SCIPreleaseCons(getSCIP(), &cons) ); + + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &sumexpr) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &terms[1]) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &terms[0]) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &xexpr) ); +} + +void ScipModelAPI::AddConstraint(const PowConstraint &cc) { + SCIP_VAR* x = getPROBDATA()->vars[cc.GetArguments()[0]]; + SCIP_VAR* res = getPROBDATA()->vars[cc.GetResultVar()]; + + assert(x != NULL); + assert(res != NULL); + + SCIP_EXPR* terms[2]; + SCIP_EXPR* xexpr; + SCIP_EXPR* sumexpr; + SCIP_Real coefs[2] = {1.0, -1.0}; + + SCIP_CCALL( SCIPcreateExprVar(getSCIP(), &xexpr, x, NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprVar(getSCIP(), &terms[0], res, NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprPow(getSCIP(), &terms[1], xexpr, cc.GetParameters()[0], NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprSum(getSCIP(), &sumexpr, 2, terms, coefs, 0.0, NULL, NULL) ); + + SCIP_CONS* cons; + SCIP_CCALL( SCIPcreateConsBasicNonlinear(getSCIP(), &cons, cc.GetName(), sumexpr, 0.0, 0.0) ); + SCIP_CCALL( SCIPaddCons(getSCIP(), cons) ); + SCIP_CCALL( SCIPreleaseCons(getSCIP(), &cons) ); + + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &sumexpr) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &terms[1]) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &terms[0]) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &xexpr) ); +} + +void ScipModelAPI::AddConstraint(const SinConstraint &cc) { + SCIP_VAR* x = getPROBDATA()->vars[cc.GetArguments()[0]]; + SCIP_VAR* res = getPROBDATA()->vars[cc.GetResultVar()]; + + assert(x != NULL); + assert(res != NULL); + + SCIP_EXPR* terms[2]; + SCIP_EXPR* xexpr; + SCIP_EXPR* sumexpr; + SCIP_Real coefs[2] = {1.0, -1.0}; + + SCIP_CCALL( SCIPcreateExprVar(getSCIP(), &xexpr, x, NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprVar(getSCIP(), &terms[0], res, NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprSin(getSCIP(), &terms[1], xexpr, NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprSum(getSCIP(), &sumexpr, 2, terms, coefs, 0.0, NULL, NULL) ); + + SCIP_CONS* cons; + SCIP_CCALL( SCIPcreateConsBasicNonlinear(getSCIP(), &cons, cc.GetName(), sumexpr, 0.0, 0.0) ); + SCIP_CCALL( SCIPaddCons(getSCIP(), cons) ); + SCIP_CCALL( SCIPreleaseCons(getSCIP(), &cons) ); + + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &sumexpr) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &terms[1]) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &terms[0]) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &xexpr) ); +} + +void ScipModelAPI::AddConstraint(const CosConstraint &cc) { + SCIP_VAR* x = getPROBDATA()->vars[cc.GetArguments()[0]]; + SCIP_VAR* res = getPROBDATA()->vars[cc.GetResultVar()]; + + assert(x != NULL); + assert(res != NULL); + + SCIP_EXPR* terms[2]; + SCIP_EXPR* xexpr; + SCIP_EXPR* sumexpr; + SCIP_Real coefs[2] = {1.0, -1.0}; + + SCIP_CCALL( SCIPcreateExprVar(getSCIP(), &xexpr, x, NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprVar(getSCIP(), &terms[0], res, NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprCos(getSCIP(), &terms[1], xexpr, NULL, NULL) ); + SCIP_CCALL( SCIPcreateExprSum(getSCIP(), &sumexpr, 2, terms, coefs, 0.0, NULL, NULL) ); + + SCIP_CONS* cons; + SCIP_CCALL( SCIPcreateConsBasicNonlinear(getSCIP(), &cons, cc.GetName(), sumexpr, 0.0, 0.0) ); + SCIP_CCALL( SCIPaddCons(getSCIP(), cons) ); + SCIP_CCALL( SCIPreleaseCons(getSCIP(), &cons) ); + + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &sumexpr) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &terms[1]) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &terms[0]) ); + SCIP_CCALL( SCIPreleaseExpr(getSCIP(), &xexpr) ); +} + +void ScipModelAPI::FinishProblemModificationPhase() { +} + + +} // namespace mp diff --git a/solvers/scip/scipmodelapi.h b/solvers/scip/scipmodelapi.h new file mode 100644 index 000000000..3ae4293d6 --- /dev/null +++ b/solvers/scip/scipmodelapi.h @@ -0,0 +1,140 @@ +#ifndef SCIPMODELAPI_H +#define SCIPMODELAPI_H + +#include + +#include "mp/env.h" +#include "scipcommon.h" +#include "mp/flat/model_api_base.h" +#include "mp/flat/constr_std.h" + +namespace mp { + +class ScipModelAPI : + public ScipCommon, public EnvKeeper, + public BasicFlatModelAPI +{ + using BaseModelAPI = BasicFlatModelAPI; + +private: + void linearHelper(const int* pvars, const double* pcoefs, const size_t size, const char* name, const double lb, const double ub); + void helpIndicatorLin(const int* pvars, const double* pcoefs, const size_t size, const char* name, const double rhs, const int binary_var, const int binary_value, SCIP_Bool lessthanineq); + void helpQuad(const char* name, const int* lt_pvars, const double* lt_pcoefs, const size_t lt_size, const int* qt_pvars1, const int* qt_pvars2, const double* qt_pcoefs, const size_t qt_size, const double lb, const double ub); + +public: + /// Construct + ScipModelAPI(Env& e) : EnvKeeper(e) { } + + /// Class name + static const char* GetTypeName() { return "ScipModelAPI"; } + + /// If any driver options added from here + void InitCustomOptions() { } + + /// Called before problem input. + /// Model info can be used to preallocate memory. + void InitProblemModificationPhase(const FlatModelInfo*); + /// After + void FinishProblemModificationPhase(); + + /// Implement the following functions using the solver's API + void AddVariables(const VarArrayDef& ); + void SetLinearObjective( int iobj, const LinearObjective& lo ); + /// Whether accepting quadratic objectives: + /// 0 - no, 1 - convex, 2 - nonconvex + static int AcceptsQuadObj() { return 0; } + void SetQuadraticObjective(int iobj, const QuadraticObjective& qo); + + //////////////////////////// GENERAL CONSTRAINTS //////////////////////////// + USE_BASE_CONSTRAINT_HANDLERS(BaseModelAPI) + + /// For each suppoted constraint type, add the ACCEPT_CONSTRAINT macro + /// and the relative AddConstraint function. + /// Below some typical constraint handlers of a MIP solver. + /// Further constraint types which could be handled natively by some solvers: + /// - IndicatorConstraint(Lin/Quad)(LE/EQ/GE) + /// - Multidirectional indicators Cond(Lin/Quad)Con(LT/LE/EQ/GE/GT), where + /// the implication direction () depends in the context + /// - Complementarity + /// - Logical, counting, piecewise-linear constraints. + /// See \a constr_std.h and other drivers. + + + /// The linear range constraint, if fully supported with basis info etc. + ACCEPT_CONSTRAINT(LinConRange, Recommended, CG_Linear) + void AddConstraint(const LinConRange& lc); + + /// LinCon(LE/EQ/GE) should have 'Recommended' for all backends + /// and have an implementation, + /// or a conversion rule is needed in a derived FlatConverter + ACCEPT_CONSTRAINT(LinConLE, Recommended, CG_Linear) + void AddConstraint(const LinConLE& lc); + ACCEPT_CONSTRAINT(LinConEQ, Recommended, CG_Linear) + void AddConstraint(const LinConEQ& lc); + ACCEPT_CONSTRAINT(LinConGE, Recommended, CG_Linear) + void AddConstraint(const LinConGE& lc); + + /// Ask if the solver accepts non-convex quadratic constraints + static constexpr bool AcceptsNonconvexQC() { return true; } + + /// QuadConRange is optional. + ACCEPT_CONSTRAINT(QuadConRange, Recommended, CG_Quadratic) + void AddConstraint(const QuadConRange& qc); + + /// If using quadratics, + /// QuadCon(LE/EQ/GE) should have 'Recommended' + /// and have an implementation. + ACCEPT_CONSTRAINT(QuadConLE, Recommended, CG_Quadratic) + void AddConstraint(const QuadConLE& qc); + ACCEPT_CONSTRAINT(QuadConEQ, Recommended, CG_Quadratic) + void AddConstraint(const QuadConEQ& qc); + ACCEPT_CONSTRAINT(QuadConGE, Recommended, CG_Quadratic) + void AddConstraint(const QuadConGE& qc); + + /// Linear indicator constraints can be used as + /// auxiliary constraints for logical conditions. + /// If not handled, the compared expressions need + /// deducible finite bounds for a big-M redefinition. + ACCEPT_CONSTRAINT(AbsConstraint, Recommended, CG_General) + void AddConstraint(const AbsConstraint& absc); + ACCEPT_CONSTRAINT(AndConstraint, AcceptedButNotRecommended, CG_General) + void AddConstraint(const AndConstraint& cc); + ACCEPT_CONSTRAINT(OrConstraint, Recommended, CG_General) + void AddConstraint(const OrConstraint& dc); + ACCEPT_CONSTRAINT(IndicatorConstraintLinLE, AcceptedButNotRecommended, CG_General) + void AddConstraint(const IndicatorConstraintLinLE& mc); + ACCEPT_CONSTRAINT(IndicatorConstraintLinEQ, AcceptedButNotRecommended, CG_General) + void AddConstraint(const IndicatorConstraintLinEQ& mc); + ACCEPT_CONSTRAINT(IndicatorConstraintLinGE, AcceptedButNotRecommended, CG_General) + void AddConstraint(const IndicatorConstraintLinGE& mc); + + /// Cones + /* + ACCEPT_CONSTRAINT(QuadraticConeConstraint, Recommended, CG_Conic) + void AddConstraint(const QuadraticConeConstraint& qc);*/ + + /// SOS constraints can be used by AMPL for redefinition of + /// piecewise-linear expressions. + /// Set ``option pl_linearize 0;`` in AMPL if the solver + /// supports PL natively. + ACCEPT_CONSTRAINT(SOS1Constraint, Recommended, CG_SOS) + void AddConstraint(const SOS1Constraint& cc); + ACCEPT_CONSTRAINT(SOS2Constraint, AcceptedButNotRecommended, CG_SOS) + void AddConstraint(const SOS2Constraint& cc); + + /// SCIP nonlinear generals + ACCEPT_CONSTRAINT(ExpConstraint, Recommended, CG_General) + void AddConstraint(const ExpConstraint& cc); + ACCEPT_CONSTRAINT(LogConstraint, Recommended, CG_General) + void AddConstraint(const LogConstraint& cc); + ACCEPT_CONSTRAINT(PowConstraint, Recommended, CG_General) + void AddConstraint(const PowConstraint& cc); + ACCEPT_CONSTRAINT(SinConstraint, Recommended, CG_General) + void AddConstraint(const SinConstraint& cc); + ACCEPT_CONSTRAINT(CosConstraint, AcceptedButNotRecommended, CG_General) //pretty slow + void AddConstraint(const CosConstraint& cc); +}; + +} // namespace mp + +#endif // SCIPMODELAPI_H diff --git a/test/end2end/cases/categorized/fast/cp_global_constraints/modellist.json b/test/end2end/cases/categorized/fast/cp_global_constraints/modellist.json index 21cf20d78..5bc18026a 100644 --- a/test/end2end/cases/categorized/fast/cp_global_constraints/modellist.json +++ b/test/end2end/cases/categorized/fast/cp_global_constraints/modellist.json @@ -70,6 +70,7 @@ "name" : "x-multmip1 [if, >]", "files": ["x-multmip1.mod", "multmip1.dat"], "tags" : ["logical"], + "options": { "solution_round": "6" }, "objective" : 218125 }, { diff --git a/test/end2end/cases/categorized/fast/general_constraints/modellist.json b/test/end2end/cases/categorized/fast/general_constraints/modellist.json index aece27d45..39042dc4e 100644 --- a/test/end2end/cases/categorized/fast/general_constraints/modellist.json +++ b/test/end2end/cases/categorized/fast/general_constraints/modellist.json @@ -22,7 +22,9 @@ { "name" : "poly_01", "tags" : ["polynomial"], - "options": { "ANYSOLVER_options": "nonconvex=2" }, + "options": { + "gurobi_options": "nonconvex=2" + }, "values": { "x": 0.1663127474012448 } @@ -30,7 +32,9 @@ { "name" : "and_or_max_poly_01", "tags" : ["logical", "polynomial"], - "options": { "ANYSOLVER_options": "nonconvex=2" }, + "options": { + "gurobi_options": "nonconvex=2" + }, "objective" : -40 } ] diff --git a/test/end2end/cases/categorized/fast/qp/modellist.json b/test/end2end/cases/categorized/fast/qp/modellist.json index 13b2218a4..98cdac256 100644 --- a/test/end2end/cases/categorized/fast/qp/modellist.json +++ b/test/end2end/cases/categorized/fast/qp/modellist.json @@ -55,7 +55,7 @@ { "name" : "ellipse_nonconvex testMIPstart", "objective" : 0.8, - "tags" : ["quadraticnonconvex"], + "tags" : ["quadraticnonconvex", "mipstart"], "options": { "gurobi_options": "nonconvex=2 debug=1" }, @@ -67,7 +67,7 @@ { "name" : "ellipse_nonconvex_abs testMIPstart", "objective" : 3.6, - "tags" : ["quadraticnonconvex"], + "tags" : ["quadraticnonconvex", "mipstart"], "options": { "gurobi_options": "nonconvex=2 debug=1" }, diff --git a/test/end2end/cases/categorized/fast/suf_common/modellist.json b/test/end2end/cases/categorized/fast/suf_common/modellist.json index a757b6ba1..3157006cf 100644 --- a/test/end2end/cases/categorized/fast/suf_common/modellist.json +++ b/test/end2end/cases/categorized/fast/suf_common/modellist.json @@ -107,6 +107,7 @@ { "name": "sos_01", "tags": [ "linear", "integer", "sos"], + "options": { "scip_options": "acc:sos2=2" }, "objective": 32, "values": { "x[2]": 3, diff --git a/test/end2end/cases/categorized/moderate/linear/modellist.json b/test/end2end/cases/categorized/moderate/linear/modellist.json index e7ba48034..50f0ac6c7 100644 --- a/test/end2end/cases/categorized/moderate/linear/modellist.json +++ b/test/end2end/cases/categorized/moderate/linear/modellist.json @@ -8,13 +8,21 @@ "name" : "cut", "objective" : 47, "tags" : ["linear", "integer", "script"], - "files" : ["cut.run"] + "files" : ["cut.run"], + "options": { + "scip_options": "pre:maxrounds=0 pro:maxrounds=0 pro:maxroundsroot=0", + "comment": "Disable presolve and propagation for SCIP for duals." + } }, { "name" : "cut2", "objective" : 47, "tags" : ["linear", "integer", "script"], - "files" : ["cut2.run"] + "files" : ["cut2.run"], + "options": { + "scip_options": "pre:maxrounds=0 pro:maxrounds=0 pro:maxroundsroot=0", + "comment": "Disable presolve and propagation for SCIP for duals." + } }, { "name" : "diet2 iis", diff --git a/test/end2end/cases/uncategorized/complementarity/modellist.json b/test/end2end/cases/uncategorized/complementarity/modellist.json index a0ac0bd63..b55622aa6 100644 --- a/test/end2end/cases/uncategorized/complementarity/modellist.json +++ b/test/end2end/cases/uncategorized/complementarity/modellist.json @@ -13,10 +13,6 @@ "name" : "econ", "tags" : ["linear", "continuous", "complementarity"], "files" : ["econ.mod", "econ.dat"] - }, { - "name" : "econ2", - "tags" : ["linear", "continuous", "complementarity"], - "files" : ["econ2.mod", "econ2.dat"] }, { "name" : "econnl", "tags" : ["nonlinear", "continuous", "complementarity"], diff --git a/test/end2end/scripts/python/Solver.py b/test/end2end/scripts/python/Solver.py index 3bf39913c..e0be9ab4c 100644 --- a/test/end2end/scripts/python/Solver.py +++ b/test/end2end/scripts/python/Solver.py @@ -733,6 +733,29 @@ def __init__(self, exeName, timeout=None, nthreads=None, otherOptions=None): ModelTags.quadratic_obj, ModelTags.sos} super().__init__(exeName, timeout, nthreads, otherOptions, stags) +class SCIPSolver(MPDirectSolver): + #def _setLPMethod(self, method : str): + # m = "simplex" if method == "SIMPLEX" else "ipm" + # return f"alg:method {m}" + + def _getAMPLOptionsName(self): + return "scip" + + def _setNThreads(self, threads): + return "" + + def __init__(self, exeName, timeout=None, nthreads=None, otherOptions=None): + stags = {ModelTags.continuous, ModelTags.integer, ModelTags.binary, + ModelTags.quadratic, + ModelTags.quadraticnonconvex, + ModelTags.sos, + ModelTags.nonlinear, + ModelTags.log, + ModelTags.trigonometric, + ModelTags.return_mipgap + } + super().__init__(exeName, timeout, nthreads, otherOptions, stags) + class HighsSolver(MPDirectSolver): def _setLPMethod(self, method : str): diff --git a/test/end2end/scripts/python/SolverCollection.py b/test/end2end/scripts/python/SolverCollection.py index 4168a6043..cda4dfbe0 100644 --- a/test/end2end/scripts/python/SolverCollection.py +++ b/test/end2end/scripts/python/SolverCollection.py @@ -41,6 +41,7 @@ def addStdSolvers(solvers: SolverCollection, binPath=""): solvers.addSolver(Solver.XPRESSDirectSolver(path.join(binPath,"xpress"))) solvers.addSolver(Solver.MosekSolver(path.join(binPath,"mosek"))) solvers.addSolver(Solver.CbcMPSolver(path.join(binPath, "cbc"))) + solvers.addSolver(Solver.SCIPSolver(path.join(binPath, "scip"))) solvers.addSolver(Solver.CPLEXODHSolver(path.join(binPath, "cplexodh"))) solvers.addSolver(Solver.GUROBIODHSolver(path.join(binPath, "gurobiodh")))