Repository: ampl/mp (language: python)

Issue:
Draft: GCG driver



Discussion / hints:

Fix patch:
diff --git a/solvers/CMakeLists.txt b/solvers/CMakeLists.txt
index 6b3bd2a6e..7fd97aa7d 100644
--- a/solvers/CMakeLists.txt
+++ b/solvers/CMakeLists.txt
@@ -179,6 +179,9 @@ function(add_ampl_solver name)
   endif()
   add_executable(${exename} ${dir}main.cc)
   target_link_libraries(${exename} PRIVATE ${static_lib})
+  if(add_ampl_solver_ADDITIONAL_DEPS)
+    add_dependencies(${exename} ${add_ampl_solver_ADDITIONAL_DEPS})
+  endif()
   install(TARGETS ${exename} DESTINATION bin)
   if(add_ampl_solver_ADDITIONAL_DEPS)
     add_dependencies(${exename} ${add_ampl_solver_ADDITIONAL_DEPS})
@@ -255,11 +258,17 @@ add_ampl_backend(cbcmp DLL_RUNTIME SHARED_LIB MODULE CBC
   LIBRARIES ${CBC_LIBS} ${CMAKE_DL_LIBS})
 
 
+add_ampl_backend(gcgmp DLL_RUNTIME SHARED_LIB MODULE GCG 
+  LIBRARIES ${GCG_LIBS} ${CMAKE_DL_LIBS}
+  ADDITIONAL_DEPS gcg)
+
 add_ampl_backend(scipmp DLL_RUNTIME SHARED_LIB MODULE SCIP 
-  LIBRARIES ${SCIP_LIBS} ${CMAKE_DL_LIBS})
+  LIBRARIES ${SCIP_LIBS} ${CMAKE_DL_LIBS}
+  ADDITIONAL_DEPS scip)
 
 add_ampl_backend(highsmp DLL_RUNTIME SHARED_LIB MODULE HIGHS 
-  LIBRARIES ${HIGHS_LIBS} ${CMAKE_DL_LIBS})
+  LIBRARIES ${HIGHS_LIBS} ${CMAKE_DL_LIBS}
+  ADDITIONAL_DEPS highs)
 
 add_ampl_backend(ortoolsmp DLL_RUNTIME SHARED_LIB MODULE ortoolsmp 
   LIBRARIES ${ortoolsmp_LIBS} ${CMAKE_DL_LIBS})
diff --git a/solvers/gcgmp/CHANGES.gcgmp.md b/solvers/gcgmp/CHANGES.gcgmp.md
new file mode 100644
index 000000000..bc8d7bbb3
--- /dev/null
+++ b/solvers/gcgmp/CHANGES.gcgmp.md
@@ -0,0 +1,5 @@
+Summary of recent updates to GCG for AMPL
+=============================================
+
+## 20220420
+- First release of mock driver
\ No newline at end of file
diff --git a/solvers/gcgmp/README.gcgmp.txt b/solvers/gcgmp/README.gcgmp.txt
new file mode 100644
index 000000000..062beb9c5
--- /dev/null
+++ b/solvers/gcgmp/README.gcgmp.txt
@@ -0,0 +1,5 @@
+gcg tutorial driver
+=======================
+
+The solver "gcg" is a mock solver, provided as an example - and as a template -
+on how to getting started using mp with FlatAPI to develop an mp-based solver interface.
\ No newline at end of file
diff --git a/solvers/gcgmp/gcgmp-ampls-c-api.h b/solvers/gcgmp/gcgmp-ampls-c-api.h
new file mode 100644
index 000000000..1aded1de5
--- /dev/null
+++ b/solvers/gcgmp/gcgmp-ampls-c-api.h
@@ -0,0 +1,33 @@
+#ifndef GCGAMPLSCAPI_H
+#define GCGAMPLSCAPI_H
+/*
+ * C API for MP/Gcg
+ */
+
+//#include "gcg.h"
+
+#include "mp/ampls-c-api.h"
+
+/*
+ * Below are Gcg-specific AMPLS API functions.
+ * They complement the 'public' AMPLS API defined in ampls-c-api.h.
+ */
+
+/// Initialize AMPLS gcg.
+
+/// @param slv_opt: a string of solver options
+/// (normally provided in the <solver>_options string).
+/// Can be NULL.
+/// @return pointer to struct AMPLS_MP_Solver to be populated.
+void*  AMPLSOpenGcg(const char* slv_opt, CCallbacks cb);
+
+/// Shut down solver instance
+void AMPLSCloseGcg(AMPLS_MP_Solver* slv);
+
+/// Extract the Gcg model handle
+void* GetGcgmodel(AMPLS_MP_Solver* slv);
+
+
+#endif // GCGAMPLSCAPI_H
+
+
diff --git a/solvers/gcgmp/gcgmp-lib.c b/solvers/gcgmp/gcgmp-lib.c
new file mode 100644
index 000000000..dc6ee5065
--- /dev/null
+++ b/solvers/gcgmp/gcgmp-lib.c
@@ -0,0 +1,30 @@
+#include "gcgmp/gcgmp-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 = AMPLSOpenGcg(slv_opt, cb);
+  if (!slv)
+    return NULL;
+  AMPLSLoadNLModel(slv, nl_filename);
+  return slv;
+}
+APIEXPORT void* AMPLgetGcgmodel(void* slv) {
+  return GetGcgmodel(slv);
+}
+
+APIEXPORT void AMPLwritesolution(AMPLS_MP_Solver* slv, const char* solFileName) {
+  AMPLSReportResults(slv, solFileName);
+}
+
+APIEXPORT void AMPLclosesolver(AMPLS_MP_Solver* slv) {
+  AMPLSCloseGcg(slv);
+}
\ No newline at end of file
diff --git a/solvers/gcgmp/gcgmp-modelapi-connect.cc b/solvers/gcgmp/gcgmp-modelapi-connect.cc
new file mode 100644
index 000000000..2c97139fe
--- /dev/null
+++ b/solvers/gcgmp/gcgmp-modelapi-connect.cc
@@ -0,0 +1,18 @@
+#include "mp/flat/redef/MIP/converter_mip.h"
+#include "mp/flat/model_api_connect.h"
+
+#include "gcgmpmodelapi.h"
+
+
+namespace mp {
+
+/// Defining the function in ...-modelapi-connect.cc
+/// for recompilation speed
+std::unique_ptr<BasicModelManager>
+CreateGcgModelMgr(GcgCommon& cc, Env& e,
+										 pre::BasicValuePresolver*& pPre) {
+	return CreateModelMgrWithFlatConverter<
+			GcgModelAPI, MIPFlatConverter >(cc, e, pPre);
+}
+
+}  // namespace mp
diff --git a/solvers/gcgmp/gcgmpbackend.cc b/solvers/gcgmp/gcgmpbackend.cc
new file mode 100644
index 000000000..d795b4f89
--- /dev/null
+++ b/solvers/gcgmp/gcgmpbackend.cc
@@ -0,0 +1,1092 @@
+#include <vector>
+#include <climits>
+#include <cfloat>
+
+#include "mp/env.h"
+#include "mp/flat/model_api_base.h"
+#include "gcgmpbackend.h"
+
+#include "gcg/class_partialdecomp.h"
+#include "gcg/class_detprobdata.h"
+#include "gcg/cons_decomp.hpp"
+
+extern "C" {
+  #include "gcgmp-ampls-c-api.h"    // Gcg AMPLS C API
+}
+#include "mp/ampls-cpp-api.h"
+
+namespace {
+
+
+bool InterruptGcg(void* prob) {
+  SCIPinterruptSolve(static_cast<SCIP*>(prob));
+  return true;
+}
+
+}  // namespace {}
+
+std::unique_ptr<mp::BasicBackend> CreateGcgBackend() {
+  return std::unique_ptr<mp::BasicBackend>{new mp::GcgBackend()};
+}
+
+
+namespace mp {
+
+/// Create Gcg Model Manager
+/// @param gc: the Gcg common handle
+/// @param e: environment
+/// @param pre: presolver to be returned,
+/// need it to convert solution data
+/// @return GcgModelMgr
+std::unique_ptr<BasicModelManager>
+CreateGcgModelMgr(GcgCommon&, Env&, pre::BasicValuePresolver*&);
+
+
+GcgBackend::GcgBackend() {
+  OpenSolver();
+
+  /// Create a ModelManager
+  pre::BasicValuePresolver* pPre;
+  auto data = CreateGcgModelMgr(*this, *this, pPre);
+  SetMM( std::move( data ) );
+  SetValuePresolver(pPre);
+
+  /// Copy env/lp to ModelAPI
+  copy_common_info_to_other();
+}
+
+GcgBackend::~GcgBackend() {
+  CloseSolver();
+}
+
+
+const char* GcgBackend::GetBackendName()
+  { return "GGCGBackend"; }
+
+std::string GcgBackend::GetSolverVersion() {
+  return fmt::format("{}.{}.{}", GCGmajorVersion(), 
+    GCGminorVersion(), GCGtechVersion());
+}
+
+
+bool GcgBackend::IsMIP() const {
+  return SCIPgetNOrigIntVars(getSCIP()) > 0 || SCIPgetNOrigBinVars(getSCIP()) > 0;
+}
+
+bool GcgBackend::IsQCP() const {
+  return false;
+}
+
+ArrayRef<double> GcgBackend::PrimalSolution() {
+  SCIP* scip = getSCIP();
+  int num_vars = NumVars();
+  std::vector<double> x(num_vars);
+  for (int i = 0; i < num_vars; i++)
+    x[i] = SCIPgetSolVal(scip, SCIPgetBestSol(scip), getPROBDATA()->vars[i]);
+  return x;
+}
+
+pre::ValueMapDbl GcgBackend::DualSolution() {
+  return {{ { CG_Linear, DualSolution_LP() } }};
+}
+
+ArrayRef<double> GcgBackend::DualSolution_LP() {
+  int num_cons = NumLinCons();
+  std::vector<double> pi(num_cons);
+ // int error = GCG_GetLpSolution(lp(), NULL, NULL, pi.data(), NULL);
+  int error = 0;
+  if (error)
+    pi.clear();
+  return pi;
+}
+
+double GcgBackend::ObjectiveValue() const {
+  return SCIPgetPrimalbound(getSCIP());
+}
+
+double GcgBackend::NodeCount() const {
+  return SCIPgetNNodes(getSCIP());
+}
+
+double GcgBackend::SimplexIterations() const {
+  return SCIPgetNPrimalLPIterations(getSCIP()) + SCIPgetNDualLPIterations(getSCIP());
+}
+
+int GcgBackend::BarrierIterations() const {
+  return SCIPgetNBarrierLPIterations(getSCIP());
+}
+
+void GcgBackend::ExportModel(const std::string &file) {
+  GCG_CCALL( SCIPwriteOrigProblem(getSCIP(), file.data(), NULL, FALSE) );
+}
+
+
+void GcgBackend::SetInterrupter(mp::Interrupter *inter) {
+  inter->SetHandler(InterruptGcg, getSCIP());
+  // TODO Check interrupter
+  //GCG_CCALL( CPXsetterminate (env(), &terminate_flag) );
+}
+
+void GcgBackend::Solve() {
+  if (!storedOptions_.exportFile_.empty())
+    ExportModel(storedOptions_.exportFile_);
+  if (!storedOptions_.paramRead_.empty())
+    GCG_CCALL( SCIPreadParams(getSCIP(), storedOptions_.paramRead_.c_str()) );
+  if (storedOptions_.heuristics_ != 0)
+    GCG_CCALL( SCIPsetHeuristics(getSCIP(), (SCIP_PARAMSETTING)storedOptions_.heuristics_, TRUE) );
+  if (storedOptions_.cuts_ != 0)
+    GCG_CCALL( SCIPsetSeparating(getSCIP(), (SCIP_PARAMSETTING)storedOptions_.cuts_, TRUE) );
+  if (storedOptions_.presolvings_ != 0)
+    GCG_CCALL( SCIPsetSeparating(getSCIP(), (SCIP_PARAMSETTING)storedOptions_.presolvings_, TRUE) );
+
+  InputDecomposition();
+  
+  GCG_CCALL(GCGsolve(getSCIP()));
+
+  WindupGCGSolve();
+}
+
+void GcgBackend::WindupGCGSolve() { }
+
+void GcgBackend::ReportResults() {
+  ReportGCGResults();
+  BaseBackend::ReportResults();
+}
+
+void GcgBackend::ReportGCGResults() {
+  SetStatus( ConvertGCGStatus() );
+  AddGCGMessages();
+  if (need_multiple_solutions())
+    ReportGCGPool();
+}
+std::vector<double> GcgBackend::getPoolSolution(int i)
+{
+  std::vector<double> vars(NumVars());
+ // GCG_CCALL(GCG_GetPoolSolution(lp(), i, NumVars(), NULL, vars.data()));
+  return vars;
+}
+double GcgBackend::getPoolObjective(int i)
+{
+  double obj;
+ // GCG_CCALL(GCG_GetPoolObjVal(lp(), i, &obj));
+  return obj;
+}
+void GcgBackend::ReportGCGPool() {
+  if (!IsMIP())
+    return;
+  int iPoolSolution = -1;
+  int nsolutions;
+  /*
+  while (++iPoolSolution < getIntAttr(GCG_INTATTR_POOLSOLS)) {
+    ReportIntermediateSolution(
+      { getPoolSolution(iPoolSolution),
+        {}, { getPoolObjective(iPoolSolution) } });
+  }
+  */
+}
+
+
+void GcgBackend::AddGCGMessages() {
+  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<int, std::string> GcgBackend::ConvertGCGStatus() {
+  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 GcgBackend::FinishOptionParsing() {
+  if (storedOptions_.outlev_ == 1)
+     SetSolverOption("display/verblevel", 4);
+  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 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}
+};
+
+static const mp::OptionValueInfo scoretype[] = {
+  { "maxwhi", "max white", 0},
+  { "border", "border area", 1},
+  { "classi", "classic", 2},
+  { "forswh", "max foreseeing white", 3},
+  { "spfwh", "ppc-max-white (default)", 4},
+  { "fawh", "max foreseeing white with aggregation info", 5},
+  { "spfawh", "ppc-max-white with aggregation info", 6},
+  { "bender", "experimental benders score", 7},
+  { "strong", "strong decomposition score", 8}
+};
+
+static const mp::OptionValueInfo mode[] = {
+  { "0", "Dantzig-Wolfe (default)", 0},
+  { "1", "Benders' decomposition", 1},
+  { "2", "No decomposition", 2}
+};
+
+void GcgBackend::InitCustomOptions() {
+
+  set_option_header(
+      "GCG Optimizer Options for AMPL\n"
+      "--------------------------------------------\n"
+      "\n"
+      "To set these options, assign a string specifying their values to the "
+      "AMPL option ``gcg_options``. For example::\n"
+      "\n"
+      "  ampl: option gcg_options 'mipgap=1e-6';\n");
+
+  AddStoredOption("tech:outlev outlev",
+    "0*/1: Whether to write GCG log lines (chatter) to stdout and to file.",
+    storedOptions_.outlev_);
+
+  AddSolverOption("tech:outlev-native outlev-native",
+    "0*/1/2/3/4/5: Whether to write GCG log lines (chatter) to stdout and to file (native output level of SCIP).",
+    "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 GCG 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");
+
+  /////////////////////// 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: 2000)",
+    "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: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: 0)",
+    "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);
+
+  ////////////////////// 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_.presolvings_, 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);
+
+
+  ////////////////////// DETECTION ///////////////////////
+  AddSolverOption("det:enabled",
+    "0/1: whether detection should be enabled? "
+    "\n"
+    "  | 0 - Detection not enabled\n"
+    "  | 1 - Detection enabled (default).",
+    "detection/enabled", 0, 1);
+
+  AddSolverOption("det:postprocess postprocess",
+    "0/1: whether postprocessing of complete decompositions should be enabled? "
+    "\n"
+    "  | 0 - Postprocessing of complete decompositions not enabled\n"
+    "  | 1 - Postprocessing of complete decompositions enabled (default).",
+    "detection/enabled", 0, 1);
+
+  AddSolverOption("det:maxrounds",
+    "Maximum number of detection loop rounds (default: 1) ",
+    "detection/maxrounds", 0, INT_MAX);
+
+  AddSolverOption("det:maxtime",
+    "Maximum detection time in seconds (default: 600) ",
+    "detection/maxrounds", 0, INT_MAX);
+
+  AddSolverOption("det:scoretype scoretype",
+    "Score calculation for comparing (partial) decompositions:\n"
+    "\n.. value-table::\n", "detection/scores/selected", scoretype, 4);
+
+  AddSolverOption("det:origprob-classificationenabled origprob-classificationenabled",
+    "0/1: whether classification for the original problem should be enabled? "
+    "\n"
+    "  | 0 - Classification for the original problem not enabled\n"
+    "  | 1 - Classification for the original problem enabled (default).",
+    "detection/origprob/classificationenabled", 0, 1);
+
+  AddSolverOption("det:origprob-enabled origprob-enabled",
+    "0/1: whether detection for the original problem should be enabled? "
+    "\n"
+    "  | 0 - Detection for the original problem not enabled\n"
+    "  | 1 - Detection for the original problem enabled (default).",
+    "detection/origprob/enabled", 0, 1);
+
+  AddSolverOption("det:classification-enabled classification-enabled",
+    "0/1: whether classification should be enabled? "
+    "\n"
+    "  | 0 - Classification not enabled\n"
+    "  | 1 - Classification enabled (default).",
+    "detection/classification/enabled", 0, 1);
+
+  AddSolverOption("det:maxnclassesperpartition maxnclassesperpartition",
+    "Maximum number of classes per partition (default: 9)",
+    "detection/classification/maxnclassesperpartition", 0, INT_MAX);
+
+  AddSolverOption("det:maxnclassesperpartitionforlargeprobs maxnclassesperpartitionforlargeprobs",
+    "Maximum number of classes per partition for large problems (nconss + nvars >= 50000) (default: 5)",
+    "detection/classification/maxnclassesperpartitionforlargeprobs", 0, INT_MAX);
+  
+  AddSolverOption("det:benders-enabled benders-enabled",
+    "0/1: whether benders detection should be enabled? "
+    "\n"
+    "  | 0 - Benders detection not enabled (default)\n"
+    "  | 1 - Benders detection enabled.",
+    "detection/benders/enabled", 0, 1);
+
+  AddSolverOption("det:benders-onlybinmaster benders-onlybinmaster",
+    "0/1: whether only decomposition with only binary variables in the master are searched? "
+    "\n"
+    "  | 0 - Not only decomposition with only binary variables in the master are searched (default)\n"
+    "  | 1 - Only decomposition with only binary variables in the master are searched.",
+    "detection/benders/onlybinmaster", 0, 1);
+
+  AddSolverOption("det:benders-onlycontsubpr benders-onlycontsubpr",
+    "0/1: whether only decomposition with only continiuous variables in the subproblems are searched? "
+    "\n"
+    "  | 0 - Not only decomposition with only continiuous variables in the subproblems are searched (default)\n"
+    "  | 1 - Only decomposition with only continiuous variables in the subproblems are searched.",
+    "detection/benders/onlycontsubpr", 0, 1);
+
+  //////////////////// RELAXING-GCG //////////////////////
+  AddSolverOption("gcg:bliss-enabled bliss-enabled",
+    "0/1: whether bliss should be used to check for identical blocks? "
+    "\n"
+    "  | 0 - Bliss should not be used to check for identical blocks\n"
+    "  | 1 - Bliss should be used to check for identical blocks (default).",
+    "relaxing/gcg/bliss/enabled", 0, 1);
+
+  AddSolverOption("gcg:aggregation",
+    "0/1: whether identical blocks should be aggregated (only for discretization approach)? "
+    "\n"
+    "  | 0 - Identical blocks should not be aggregated (only for discretization approach)\n"
+    "  | 1 - Identical blocks should be aggregated (only for discretization approach) (default).",
+    "relaxing/gcg/aggregation", 0, 1);
+
+  AddSolverOption("gcg:discretization discretization",
+    "0/1: whether discretization (TRUE) or convexification (FALSE) approach should be used? "
+    "\n"
+    "  | 0 - convexification approach should be used\n"
+    "  | 1 - discretization approach should be used (default).",
+    "relaxing/gcg/discretization", 0, 1);
+
+  AddSolverOption("gcg:mipdiscretization mipdiscretization",
+    "0/1: whether discretization (TRUE) or convexification (FALSE) approach should be used in mixed-integer programs?"
+    "\n"
+    "  | 0 - convexification approach should be used in mixed-integer programs\n"
+    "  | 1 - discretization approach should be used in mixed-integer programs (default).",
+    "relaxing/gcg/discretization", 0, 1);
+
+  AddSolverOption("gcg:mode mode",
+    "The decomposition mode that GCG will use:\n"
+    "\n.. value-table::\n", "relaxing/gcg/mode", mode, 0);
+}
+
+void GcgBackend::InputDecomposition() {
+  bool is_presolved = SCIPgetStage(getSCIP()) >= SCIP_STAGE_PRESOLVED;
+  GCG_CCALL( SCIPallocClearMemory(scip, &getPROBDATA()->decomp) );
+  gcg::PARTIALDECOMP* decomp = new gcg::PARTIALDECOMP(getSCIP(), !is_presolved);
+  getPROBDATA()->decomp = decomp;
+
+  if (auto block0 = ReadModelSuffixInt({"block", suf::Kind::CON_BIT | suf::Kind::VAR_BIT})) {
+    auto block = GetValuePresolver().PresolveGenericInt( block0 );
+
+    ArrayRef<int> blockconss = block.GetConValues()(CG_Linear);
+    for (size_t i = 0; i < blockconss.size(); i++) {
+      getPROBDATA()->decomp->fixConsToBlock(getPROBDATA()->linconss[i], blockconss[i]);
+    }
+
+    ArrayRef<int> blockvars = block.GetVarValues()();
+    for (size_t i = 0; i < blockvars.size(); i++) {
+      int varindex = getPROBDATA()->decomp->getDetprobdata()->getIndexForVar(getPROBDATA()->vars[i]);
+      getPROBDATA()->decomp->fixVarToBlock(varindex, blockvars[i]);
+    }
+  }
+
+  if (auto master0 = ReadModelSuffixInt({"master", suf::Kind::CON_BIT | suf::Kind::VAR_BIT})) {
+    auto master = GetValuePresolver().PresolveGenericInt( master0 );
+
+    ArrayRef<int> masterconss = master.GetConValues()(CG_Linear);
+    for (size_t i = 0; i < masterconss.size(); i++) {
+      if (masterconss[i] == 1)
+        getPROBDATA()->decomp->fixConsToMaster(getPROBDATA()->linconss[i]);
+    }
+
+    ArrayRef<int> mastervars = master.GetVarValues()();
+    for (size_t i = 0; i < mastervars.size(); i++) {
+      if (mastervars[i] == 1) {
+        int varindex = getPROBDATA()->decomp->getDetprobdata()->getIndexForVar(getPROBDATA()->vars[i]);
+        getPROBDATA()->decomp->fixVarToMaster(varindex);
+      }
+    }
+  }
+
+  if (auto linking0 = ReadModelSuffixInt({"linking", suf::Kind::VAR_BIT})) {
+    auto linking = GetValuePresolver().PresolveGenericInt( linking0 );
+
+    ArrayRef<int> linkingvars = linking.GetConValues()(CG_Linear);
+    for (size_t i = 0; i < linkingvars.size(); i++) {
+      if (linkingvars[i] == 1) {
+        int varindex = getPROBDATA()->decomp->getDetprobdata()->getIndexForVar(getPROBDATA()->vars[i]);
+        getPROBDATA()->decomp->fixVarToLinking(varindex);
+      }
+    }
+  }
+
+  if (getPROBDATA()->decomp->getNOpenconss() != getPROBDATA()->decomp->getNConss() || getPROBDATA()->decomp->getNOpenvars() != getPROBDATA()->decomp->getNVars()) {
+    getPROBDATA()->decomp->prepare();
+    if (getPROBDATA()->decomp->isComplete())
+      getPROBDATA()->decomp->setUsergiven(gcg::USERGIVEN::COMPLETE);
+    else
+      getPROBDATA()->decomp->setUsergiven(gcg::USERGIVEN::PARTIAL);
+    GCGconshdlrDecompAddPreexisitingPartialDec(getSCIP(), getPROBDATA()->decomp);
+    GCGpresolve(getSCIP());
+  }
+}
+
+double GcgBackend::MIPGap() {
+  return SCIPgetGap(getSCIP())<Infinity() ? SCIPgetGap(getSCIP()) : AMPLInf();
+}
+double GcgBackend::BestDualBound() {
+  return SCIPgetDualbound(getSCIP());
+}
+
+double GcgBackend::MIPGapAbs() {
+  double gapabs = std::fabs(ObjectiveValue() - BestDualBound());
+  return gapabs<Infinity() ? gapabs : AMPLInf();
+}
+
+
+ArrayRef<int> GcgBackend::VarStatii() {
+  
+  std::vector<int> vars(NumVars());
+  /*
+  GCG_GetBasis(lp(), vars.data(), NULL);
+  for (auto& s : vars) {
+    switch (s) {
+    case GCG_BASIS_BASIC:
+      s = (int)BasicStatus::bas;
+      break;
+    case GCG_BASIS_LOWER:
+      s = (int)BasicStatus::low;
+      break;
+    case GCG_BASIS_UPPER:
+      s = (int)BasicStatus::upp;
+      break;
+    case GCG_BASIS_SUPERBASIC:
+      s = (int)BasicStatus::sup;
+      break;
+    case GCG_BASIS_FIXED:
+      s = (int)BasicStatus::equ;
+      break;
+    default:
+      MP_RAISE(fmt::format("Unknown Gcg VBasis value: {}", s));
+    }
+  }
+  */
+  return vars;
+}
+
+ArrayRef<int> GcgBackend::ConStatii() {
+
+  std::vector<int> cons(NumLinCons());
+  /*
+  GCG_GetBasis(lp(), NULL, cons.data());
+  for (auto& s : cons) {
+    switch (s) {
+    case GCG_BASIS_BASIC:
+      s = (int)BasicStatus::bas;
+      break;
+    case GCG_BASIS_LOWER:
+      s = (int)BasicStatus::low;
+      break;
+    case GCG_BASIS_UPPER:
+      s = (int)BasicStatus::upp;
+      break;
+    case GCG_BASIS_SUPERBASIC:
+      s = (int)BasicStatus::sup;
+      break;
+    case GCG_BASIS_FIXED:
+      s = (int)BasicStatus::equ;
+      break;
+    default:
+      MP_RAISE(fmt::format("Unknown Gcg VBasis value: {}", s));
+    }
+  }*/
+  return cons;
+}
+
+void GcgBackend::VarStatii(ArrayRef<int> vst) {
+  int index[1];
+  std::vector<int> stt(vst.data(), vst.data() + vst.size());
+  /*
+  for (auto j = stt.size(); j--; ) {
+    auto& s = stt[j];
+    switch ((BasicStatus)s) {
+    case BasicStatus::bas:
+      s = GCG_BASIS_BASIC;
+      break;
+    case BasicStatus::low:
+      s = GCG_BASIS_LOWER;
+      break;
+    case BasicStatus::equ:
+      s = GCG_BASIS_FIXED;
+      break;
+    case BasicStatus::upp:
+      s = GCG_BASIS_UPPER;
+      break;
+    case BasicStatus::sup:
+    case BasicStatus::btw:
+      s = GCG_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(!GCG_GetColInfo(lp(), GCG_DBLINFO_LB, 1, index, &lb) && 
+        !GCG_GetColInfo(lp(), GCG_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));
+    }
+  }
+  GCG_SetBasis(lp(), stt.data(), NULL);
+  */
+}
+
+void GcgBackend::ConStatii(ArrayRef<int> cst) {
+  /*
+  std::vector<int> stt(cst.data(), cst.data() + cst.size());
+  for (auto& s : stt) {
+    switch ((BasicStatus)s) {
+    case BasicStatus::bas:
+      s = GCG_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 = GCG_BASIS_SUPERBASIC;
+      break;
+    default:
+      MP_RAISE(fmt::format("Unknown AMPL con status value: {}", s));
+    }
+  }
+  GCG_SetBasis(lp(), NULL, stt.data());
+  */
+}
+
+SolutionBasis GcgBackend::GetBasis() {
+  std::vector<int> varstt = VarStatii();
+  std::vector<int> 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 GcgBackend::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 GcgBackend::ComputeIIS() {
+  //GCG_CCALL(GCG_ComputeIIS(lp()));
+  SetStatus(ConvertGCGStatus());   // could be new information
+}
+
+IIS GcgBackend::GetIIS() {
+  auto variis = VarsIIS();
+  auto coniis = ConsIIS();
+  auto mv = GetValuePresolver().PostsolveIIS(
+    { variis, coniis });
+  return { mv.GetVarValues()(), mv.GetConValues()() };
+}
+
+ArrayRef<int> GcgBackend::VarsIIS() {
+  return ArrayRef<int>();
+//  return getIIS(lp(), NumVars(), GCG_GetColLowerIIS, GCG_GetColUpperIIS);
+}
+pre::ValueMapInt GcgBackend::ConsIIS() {
+  /*auto iis_lincon = getIIS(lp(), NumLinCons(), GCG_GetRowLowerIIS, GCG_GetRowUpperIIS);
+
+  std::vector<int> iis_soscon(NumSOSCons());
+  GCG_GetSOSIIS(lp(), NumSOSCons(), NULL, iis_soscon.data());
+  ConvertIIS2AMPL(iis_soscon);
+
+  std::vector<int> iis_indicon(NumIndicatorCons());
+  GCG_GetIndicatorIIS(lp(), NumIndicatorCons(), NULL, iis_indicon.data());
+  ConvertIIS2AMPL(iis_indicon);
+
+  return { {{ CG_Linear, iis_lincon },
+      { CG_SOS, iis_soscon },
+      { CG_Logical, iis_indicon }} };
+      */
+  return { {{ 0, std::vector<int>()}} };
+}
+
+void GcgBackend::AddMIPStart(ArrayRef<double> x0, ArrayRef<int> sparsity) {
+  //GCG_CCALL(GCG_AddMipStart(lp(), NumVars(), NULL, const_cast<double*>(x0.data())));
+
+}
+
+
+} // namespace mp
+
+
+// AMPLs
+void* AMPLSOpenGcg(
+  const char* slv_opt, CCallbacks cb = {}) {
+  return AMPLS__internal__Open(std::unique_ptr<mp::BasicBackend>{new mp::GcgBackend()},
+    slv_opt, cb);
+}
+
+void AMPLSCloseGcg(AMPLS_MP_Solver* slv) {
+  AMPLS__internal__Close(slv);
+}
+
+void* GetGcgmodel(AMPLS_MP_Solver* slv) {
+  return
+    dynamic_cast<mp::GcgBackend*>(AMPLSGetBackend(slv))->getSCIP();
+}
diff --git a/solvers/gcgmp/gcgmpbackend.h b/solvers/gcgmp/gcgmpbackend.h
new file mode 100644
index 000000000..4d3c82c3e
--- /dev/null
+++ b/solvers/gcgmp/gcgmpbackend.h
@@ -0,0 +1,180 @@
+#ifndef MP_GCG_BACKEND_H_
+#define MP_GCG_BACKEND_H_
+
+#include <string>
+
+#include "mp/backend-mip.h"
+#include "mp/flat/backend_flat.h"
+#include "gcgmpcommon.h"
+
+namespace mp {
+
+class GcgBackend :
+    public FlatBackend< MIPBackend<GcgBackend> >,
+    public GcgCommon
+{
+  using BaseBackend = FlatBackend< MIPBackend<GcgBackend> >;
+
+  //////////////////// [[ The public interface ]] //////////////////////
+public:
+  GcgBackend();
+  ~GcgBackend();
+
+  /// Prefix used for the <prefix>_options environment variable
+  static const char* GetAMPLSolverName() { return "gcg"; }
+
+  /// AMPL driver name displayed in messages
+  static const char* GetAMPLSolverLongName() { return "AMPL-GCG"; }
+  /// Solver name displayed in messages
+  static const char* GetSolverName() { return "GCG"; }
+  /// 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, true)
+
+  /**
+  * Get/Set AMPL var/con statii
+  **/
+  ALLOW_STD_FEATURE(BASIS, true)
+  // TODO If getting/setting a basis is supported, implement the 
+  // accessor and the setter below
+  SolutionBasis GetBasis() override;
+  void SetBasis(SolutionBasis) override;
+
+  /**
+  * MIP warm start
+  **/
+  // TODO 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<double> x0,
+									 ArrayRef<int> sparsity) override;
+
+
+ /**
+  * Get MIP Gap
+  **/
+  // TODO Implement to 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
+  **/
+  // TODO Implement to return the best dual bound value
+  // (adds option mip:bestbound)
+  ALLOW_STD_FEATURE(RETURN_BEST_DUAL_BOUND, true)
+  double BestDualBound() override;
+
+  /**
+  * Compute the IIS and obtain relevant values
+  **/
+  ALLOW_STD_FEATURE(IIS, false)
+  /// Compute IIS
+  void ComputeIIS() override;
+  /// Retrieve IIS elements
+  IIS GetIIS() 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<double> GetObjectiveValues() override
+  { return std::vector<double>{ObjectiveValue()}; } 
+
+
+  //////////////////// [[ Implementation details ]] //////////////////////
+  ///////////////////////////////////////////////////////////////////////////////
+protected:
+  void ExportModel(const std::string& file);
+
+  double ObjectiveValue() const;
+
+  /// Solution values. The vectors are emptied if not available
+  ArrayRef<double> PrimalSolution() override;
+  pre::ValueMapDbl DualSolution() override;
+  ArrayRef<double> DualSolution_LP();
+
+  void WindupGCGSolve();
+
+  void ReportResults() override;
+  void ReportGCGResults();
+
+  void ReportGCGPool();
+
+  std::vector<double> getPoolSolution(int i);
+  double getPoolObjective(int i);
+
+  /// Solution attributes
+  double NodeCount() const;
+  double SimplexIterations() const;
+  int BarrierIterations() const;
+
+  std::pair<int, std::string> ConvertGCGStatus();
+  void AddGCGMessages();
+
+  ArrayRef<int> VarStatii();
+  ArrayRef<int> ConStatii();
+  void VarStatii(ArrayRef<int>);
+  void ConStatii(ArrayRef<int>);
+
+  ArrayRef<int> VarsIIS();
+  pre::ValueMapInt ConsIIS();
+
+
+private:
+  void InputDecomposition();
+  
+  /// These options are stored in the class
+  struct Options {
+    std::string exportFile_, logFile_, paramRead_;
+    int heuristics_ = 0;
+    int cuts_ = 0;
+    int presolvings_ = 0;
+    int outlev_ = 0;
+  };
+  Options storedOptions_;
+
+
+};
+
+}  // namespace mp
+
+#endif  // MP_GCG_BACKEND_H_
diff --git a/solvers/gcgmp/gcgmpcommon.cc b/solvers/gcgmp/gcgmpcommon.cc
new file mode 100644
index 000000000..b41ee6b18
--- /dev/null
+++ b/solvers/gcgmp/gcgmpcommon.cc
@@ -0,0 +1,152 @@
+#include "mp/format.h"
+#include "gcgmpcommon.h"
+
+#include "gcg/gcgplugins.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 )
+  {
+    GCG_CCALL( SCIPreleaseVar(scip, &(*probdata)->vars[i]) );
+  }
+  SCIPfreeBlockMemoryArray(scip, &(*probdata)->vars, (*probdata)->nvars);
+
+  SCIPfreeMemory(scip, probdata);
+
+  return SCIP_OKAY;
+}
+namespace mp {
+
+void GcgCommon::OpenSolver() {
+  int status = 0;
+  SCIP* scip = NULL;
+  SCIP_PROBDATA* probdata = NULL;
+
+  // initialize SCIP
+  status = SCIPcreate(&scip);
+  setSCIP(scip); // Assign it
+
+  // include default GCG plugins
+  GCG_CCALL( SCIPincludeGcgPlugins(scip) );
+
+  // initialize empty SCIP problem
+  GCG_CCALL( SCIPallocClearMemory(scip, &probdata) );
+  GCG_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 GcgCommon::CloseSolver() {
+  SCIP* scip = getSCIP();
+
+  // free SCIP
+  GCG_CCALL( SCIPfree(&scip) );
+}
+
+int GcgCommon::NumLinCons() const {
+  return getPROBDATA()->nlinconss;
+}
+
+int GcgCommon::NumVars() const {
+  return getPROBDATA()->nvars;
+}
+
+int GcgCommon::NumObjs() const {
+  return 1;
+}
+
+int GcgCommon::NumQPCons() const {
+  return 0;
+}
+
+int GcgCommon::NumSOSCons() const {
+  return 0;
+}
+
+int GcgCommon::NumIndicatorCons() const {
+  return 0;
+}
+
+
+void GcgCommon::GetSolverOption(const char* key, int &value) const {
+  if (SCIPparamGetType(SCIPgetParam(getSCIP(), key))==SCIP_PARAMTYPE_BOOL) {
+    SCIP_Bool buffer;
+    GCG_CCALL( SCIPgetBoolParam(getSCIP(), key, &buffer) );
+    value = (int)buffer;
+  }
+  else
+    GCG_CCALL( SCIPgetIntParam(getSCIP(), key, &value) );
+}
+
+void GcgCommon::SetSolverOption(const char* key, int value) {
+  if (SCIPparamGetType(SCIPgetParam(getSCIP(), key))==SCIP_PARAMTYPE_BOOL)
+    GCG_CCALL( SCIPsetBoolParam(getSCIP(), key, value) );
+  else
+    GCG_CCALL( SCIPsetIntParam(getSCIP(), key, value) );
+}
+
+void GcgCommon::GetSolverOption(const char* key, double &value) const {
+  GCG_CCALL( SCIPgetRealParam(getSCIP(), key, &value) );
+}
+
+void GcgCommon::SetSolverOption(const char* key, double value) {
+  GCG_CCALL( SCIPsetRealParam(getSCIP(), key, value) );
+}
+
+void GcgCommon::GetSolverOption(const char* key, std::string &value) const {
+  char* buffer;
+  GCG_CCALL( SCIPallocBufferArray(getSCIP(), &buffer, SCIP_MAXSTRLEN) );
+  if (SCIPparamGetType(SCIPgetParam(getSCIP(), key))==SCIP_PARAMTYPE_CHAR)
+    GCG_CCALL( SCIPgetCharParam(getSCIP(), key, buffer) );
+  else
+    GCG_CCALL( SCIPgetStringParam(getSCIP(), key, &buffer) );
+  value = buffer;
+  SCIPfreeBufferArray(getSCIP(), &buffer);
+}
+
+void GcgCommon::SetSolverOption(const char* key, const std::string& value) {
+  if (SCIPparamGetType(SCIPgetParam(getSCIP(), key))==SCIP_PARAMTYPE_CHAR)
+    GCG_CCALL( SCIPsetCharParam(getSCIP(), key, value.c_str()[0]) );
+  else
+    GCG_CCALL( SCIPsetStringParam(getSCIP(), key, value.c_str()) );
+}
+
+
+double GcgCommon::Infinity() const { 
+  double inf;
+  GetSolverOption("numerics/infinity", inf);
+  return inf;  
+}
+
+double GcgCommon::MinusInfinity() const { 
+  double inf;
+  GetSolverOption("numerics/infinity", inf);
+  return -inf;  
+}
+
+bool GcgCommon::IsContinuous() {
+  // True iff gcg 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/gcgmp/gcgmpcommon.h b/solvers/gcgmp/gcgmpcommon.h
new file mode 100644
index 000000000..b806a4796
--- /dev/null
+++ b/solvers/gcgmp/gcgmpcommon.h
@@ -0,0 +1,87 @@
+#ifndef GCGCOMMON_H
+#define GCGCOMMON_H
+
+#include <string>
+
+#include "mp/backend-to-model-api.h"
+
+#include "scip/scip.h"
+#include "gcg/gcg.h"
+
+#include "gcg/class_partialdecomp.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 */
+
+   gcg::PARTIALDECOMP*   decomp;             /**< user partialdecomp */
+};
+
+namespace mp {
+
+/// Information shared by both
+/// `GcgBackend` and `GcgModelAPI`
+struct GcgCommonInfo {
+  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 Gcg classes
+class GcgCommon :
+    public Backend2ModelAPIConnector<GcgCommonInfo> {
+public:
+  /// These methods access Gcg 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);
+
+  /// GCG own infinity
+  double Infinity() const;
+  double MinusInfinity() const;
+
+  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;
+};
+
+
+/// Convenience macro
+// 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 GCG_RETCODE_OK 1
+#define GCG_CCALL( call ) do { if (int e = (call) != GCG_RETCODE_OK) \
+  throw std::runtime_error( \
+    fmt::format("  Call failed: '{}' with code {}", #call, e )); } while (0)
+
+} // namespace mp
+
+#endif // GCGCOMMON_H
diff --git a/solvers/gcgmp/gcgmpmodelapi.cc b/solvers/gcgmp/gcgmpmodelapi.cc
new file mode 100644
index 000000000..94f09b5cf
--- /dev/null
+++ b/solvers/gcgmp/gcgmpmodelapi.cc
@@ -0,0 +1,97 @@
+#include "gcgmpmodelapi.h"
+
+
+namespace mp {
+
+void GcgModelAPI::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;
+  GCG_CCALL( SCIPallocBlockMemoryArray(getSCIP(), &getPROBDATA()->linconss, getPROBDATA()->nlinconss) );
+}
+
+void GcgModelAPI::AddVariables(const VarArrayDef& v) {
+  assert( SCIPgetStage(getSCIP()) == SCIP_STAGE_PROBLEM );
+
+  getPROBDATA()->nvars = v.size();
+  GCG_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;
+    if (v.pnames() != NULL)
+      GCG_CCALL( SCIPcreateVarBasic(getSCIP(), &var, v.pnames()[i], lb, ub, objcoef, vartype) );
+    else
+      GCG_CCALL( SCIPcreateVarBasic(getSCIP(), &var, NULL, lb, ub, objcoef, vartype) );
+    GCG_CCALL( SCIPaddVar(getSCIP(), var) );
+    getPROBDATA()->vars[i] = var;
+  }
+}
+
+void GcgModelAPI::SetLinearObjective( int iobj, const LinearObjective& lo ) {
+  if (iobj<1) {
+    GCG_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++) {
+      GCG_CCALL( SCIPchgVarObj(getSCIP(), vars[lo.vars()[i]], lo.coefs()[i]) );
+    }
+  } 
+  else {
+    throw std::runtime_error("Multiple linear objectives not supported");
+  }
+}
+
+
+void GcgModelAPI::SetQuadraticObjective(int iobj, const QuadraticObjective& qo) {
+  if (1 > iobj) {
+    throw std::runtime_error("Quadratic objectives not supported");
+  }
+  else {
+    throw std::runtime_error("Multiple quadratic objectives not supported");
+  }
+}
+
+
+void GcgModelAPI::linearHelper(const int* pvars, const double* pcoefs, const size_t size, const char* name, const double lb, const double ub) {
+  SCIP_VAR** vars = NULL;
+  GCG_CCALL( SCIPallocBufferArray(getSCIP(), &vars, size) );
+  for (size_t i = 0; i < size; i++)
+    vars[i] = getPROBDATA()->vars[pvars[i]];
+
+  SCIP_CONS* cons;
+  GCG_CCALL( SCIPcreateConsBasicLinear(getSCIP(), &cons, name, (int)size, vars, (double*)pcoefs, lb, ub) );
+  GCG_CCALL( SCIPaddCons(getSCIP(), cons) );
+  getPROBDATA()->linconss[getPROBDATA()->i] = cons;
+  getPROBDATA()->i++;
+
+  SCIPfreeBufferArray(getSCIP(), &vars);
+}
+void GcgModelAPI::AddConstraint(const LinConRange& lc) {
+  linearHelper(lc.pvars(), lc.pcoefs(), lc.size(), lc.GetName(), lc.lb(), lc.ub());
+}
+void GcgModelAPI::AddConstraint(const LinConLE& lc) {
+  linearHelper(lc.pvars(), lc.pcoefs(), lc.size(), lc.GetName(), MinusInfinity(), lc.ub());
+}
+void GcgModelAPI::AddConstraint(const LinConEQ& lc) {
+  linearHelper(lc.pvars(), lc.pcoefs(), lc.size(), lc.GetName(), lc.ub(), lc.ub());
+}
+void GcgModelAPI::AddConstraint(const LinConGE& lc) {
+linearHelper(lc.pvars(), lc.pcoefs(), lc.size(), lc.GetName(), lc.lb(), Infinity());
+}
+
+void GcgModelAPI::FinishProblemModificationPhase() {
+}
+
+
+} // namespace mp
diff --git a/solvers/gcgmp/gcgmpmodelapi.h b/solvers/gcgmp/gcgmpmodelapi.h
new file mode 100644
index 000000000..5b4c0317f
--- /dev/null
+++ b/solvers/gcgmp/gcgmpmodelapi.h
@@ -0,0 +1,83 @@
+#ifndef GCGMODELAPI_H
+#define GCGMODELAPI_H
+
+#include <memory>
+
+#include "scip/cons_linear.h"
+
+#include "mp/env.h"
+#include "gcgmpcommon.h"
+#include "mp/flat/model_api_base.h"
+#include "mp/flat/constr_std.h"
+
+namespace mp {
+
+class GcgModelAPI :
+    public GcgCommon, 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);
+
+public:
+  /// Construct
+  GcgModelAPI(Env& e) : EnvKeeper(e) { }
+
+  /// Class name
+  static const char* GetTypeName() { return "GcgModelAPI"; }
+
+  /// 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();
+
+  /// TODO 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)
+
+  /// TODO 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 false; }
+};
+
+} // namespace mp
+
+#endif // GCGMODELAPI_H
diff --git a/solvers/gcgmp/main.cc b/solvers/gcgmp/main.cc
new file mode 100644
index 000000000..3fc2ba955
--- /dev/null
+++ b/solvers/gcgmp/main.cc
@@ -0,0 +1,9 @@
+#include "mp/backend-app.h"
+
+/// Declare a backend factory
+std::unique_ptr<mp::BasicBackend> CreateGcgBackend();
+
+extern "C" int main(int, char **argv) {
+  return
+      mp::RunBackendApp(argv, CreateGcgBackend);
+}
diff --git a/solvers/gcgmp/model-mgr-with-std-pb.cc b/solvers/gcgmp/model-mgr-with-std-pb.cc
new file mode 100644
index 000000000..2ffcce081
--- /dev/null
+++ b/solvers/gcgmp/model-mgr-with-std-pb.cc
@@ -0,0 +1,8 @@
+/**
+ * Generate ModelManagerWithPB<mp::Problem>
+ *
+ * Having a separate .cc should improve compilation speed
+ */
+
+#include "mp/model-mgr-with-std-pb.hpp"
+
diff --git a/solvers/scipmp/scipmpbackend.cc b/solvers/scipmp/scipmpbackend.cc
index 2d747225a..155283828 100644
--- a/solvers/scipmp/scipmpbackend.cc
+++ b/solvers/scipmp/scipmpbackend.cc
@@ -244,6 +244,8 @@ std::pair<int, std::string> ScipBackend::ConvertSCIPStatus() {
 
 
 void ScipBackend::FinishOptionParsing() {
+  if (storedOptions_.outlev_ == 1)
+    SetSolverOption("display/verblevel", 4);
   int v=-1;
   GetSolverOption("display/verblevel", v);
   set_verbose_mode(v>0);
@@ -312,8 +314,12 @@ void ScipBackend::InitCustomOptions() {
     "\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.",
+  AddStoredOption("tech:outlev outlev",
+    "0*/1: Whether to write SCIP log lines (chatter) to stdout and to file.",
+    storedOptions_.outlev_);
+
+  AddSolverOption("tech:outlev-native outlev-native",
+    "0*/1/2/3/4/5: Whether to write SCIP log lines (chatter) to stdout and to file (native output level of SCIP).",
     "display/verblevel", 0, 5);
 
   AddStoredOption("tech:exportfile writeprob writemodel",
diff --git a/solvers/scipmp/scipmpbackend.h b/solvers/scipmp/scipmpbackend.h
index a1c351019..b72f8a687 100644
--- a/solvers/scipmp/scipmpbackend.h
+++ b/solvers/scipmp/scipmpbackend.h
@@ -160,6 +160,7 @@ class ScipBackend :
     int heuristics_ = 0;
     int cuts_ = 0;
     int presolvings_ = 0;
+    int outlev_ = 0;
   };
   Options storedOptions_;
 
diff --git a/solvers/scipmp/scipmpmodelapi.cc b/solvers/scipmp/scipmpmodelapi.cc
index 67596f0b6..af601fe78 100644
--- a/solvers/scipmp/scipmpmodelapi.cc
+++ b/solvers/scipmp/scipmpmodelapi.cc
@@ -28,10 +28,10 @@ void ScipModelAPI::AddVariables(const VarArrayDef& v) {
       vartype = SCIP_VARTYPE_INTEGER;
     else
       vartype = SCIP_VARTYPE_CONTINUOUS;
-    const char* name = NULL;
     if (v.pnames() != NULL)
-      const char* name = v.pnames()[i];
-    SCIP_CCALL( SCIPcreateVarBasic(getSCIP(), &var, name, lb, ub, objcoef, vartype) );
+      SCIP_CCALL( SCIPcreateVarBasic(getSCIP(), &var, v.pnames()[i], lb, ub, objcoef, vartype) );
+    else
+      SCIP_CCALL( SCIPcreateVarBasic(getSCIP(), &var, NULL, lb, ub, objcoef, vartype) );
     SCIP_CCALL( SCIPaddVar(getSCIP(), var) );
     getPROBDATA()->vars[i] = var;
   }


Test patch:
diff --git a/test/end2end/cases/categorized/fast/complementarity/modellist.json b/test/end2end/cases/categorized/fast/complementarity/modellist.json
index eb78c884c..51fd3b9d6 100644
--- a/test/end2end/cases/categorized/fast/complementarity/modellist.json
+++ b/test/end2end/cases/categorized/fast/complementarity/modellist.json
@@ -3,6 +3,9 @@
     "name" : "econ2 compl_range",
     "files" : ["econ2.mod", "econ2.dat"],
     "tags" : ["linear", "continuous", "complementarity"],
+    "options": {
+      "gcg_options": "cvt:bigm=1e5"
+    },
     "values": {
       "Price['AA1']": 0.0,
       "Price['AC1']": 0.0,
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 a46195cb8..2db2012b0 100644
--- a/test/end2end/cases/categorized/fast/cp_global_constraints/modellist.json
+++ b/test/end2end/cases/categorized/fast/cp_global_constraints/modellist.json
@@ -11,7 +11,8 @@
     "tags": [ "logical" ],
     "objective": 22,
     "options": {
-      "cbc_options": "cvt:mip:bigM=100000"
+      "cbc_options": "cvt:mip:bigM=100000",
+      "gcg_options": "cvt:mip:bigM=100000 mode=2"
     }
   },
   {
@@ -20,7 +21,8 @@
     "tags": [ "logical" ],
     "objective": 22,
     "options": {
-      "cbc_options": "cvt:mip:bigM=100000"
+      "cbc_options": "cvt:mip:bigM=100000",
+      "gcg_options": "cvt:mip:bigM=100000 mode=2"
     }
   },
   {
@@ -57,7 +59,7 @@
     "name" : "x-multmip3 [count, if]",
     "files": ["x-multmip3.mod", "multmip3.dat"],
     "tags" : ["logical"],
-    "options": { "solution_round": "6" },
+    "options": { "solution_round": "6", "gcg_options": "mode=2" },
     "comment_options": "For solution_round, see #200",
     "objective" : 235625
   },
@@ -76,7 +78,10 @@
     "name" : "x-multmip1 [if, >]",
     "files": ["x-multmip1.mod", "multmip1.dat"],
     "tags" : ["logical"],
-    "options": { "solution_round": "6" },
+    "options": { 
+      "solution_round": "6",
+      "gcg_options": "cvt:mip:bigM=100000" 
+    },
     "comment_options": "For solution_round, see #200",
     "objective" : 218125
   },
diff --git a/test/end2end/cases/categorized/fast/logical/modellist.json b/test/end2end/cases/categorized/fast/logical/modellist.json
index 102baa525..d3a9dae2b 100644
--- a/test/end2end/cases/categorized/fast/logical/modellist.json
+++ b/test/end2end/cases/categorized/fast/logical/modellist.json
@@ -77,7 +77,10 @@
   {
     "name" : "impl_else_02",
     "objective" : -214,
-    "tags" : ["logical"]
+    "tags" : ["logical"],
+    "options": {
+      "gcg_options": "mode=2"
+    }
   },
   {
     "name" : "int_ne_01",
@@ -104,11 +107,17 @@
   },
   {
     "name" : "booleq_01",
+    "options": {
+      "gcg_options": "cvt:bigm=1e5"
+    },
     "objective" : 1,
     "tags" : ["logical"]
   },
   {
     "name" : "booleq_02",
+    "options": {
+      "gcg_options": "cvt:bigm=1e5"
+    },
     "objective" : 1,
     "tags" : ["logical"]
   }
diff --git a/test/end2end/cases/categorized/moderate/linear/modellist.json b/test/end2end/cases/categorized/moderate/linear/modellist.json
index 50f0ac6c7..853f8c89c 100644
--- a/test/end2end/cases/categorized/moderate/linear/modellist.json
+++ b/test/end2end/cases/categorized/moderate/linear/modellist.json
@@ -10,7 +10,8 @@
     "tags" : ["linear", "integer", "script"],
     "files" : ["cut.run"],
     "options": {
-      "scip_options": "pre:maxrounds=0 pro:maxrounds=0 pro:maxroundsroot=0",
+      "scip_options": "pre:settings=3 pro:maxrounds=0 pro:maxroundsroot=0 heu:settings=3",
+      "gcg_options": "pre:settings=3 pro:maxrounds=0 pro:maxroundsroot=0 heu:settings=3",
       "comment": "Disable presolve and propagation for SCIP for duals."
     }
   },
@@ -20,7 +21,8 @@
     "tags" : ["linear", "integer", "script"],
     "files" : ["cut2.run"],
     "options": {
-      "scip_options": "pre:maxrounds=0 pro:maxrounds=0 pro:maxroundsroot=0",
+      "scip_options": "pre:settings=3 pro:maxrounds=0 pro:maxroundsroot=0 heu:settings=3",
+      "gcg_options": "pre:settings=3 pro:maxrounds=0 pro:maxroundsroot=0 heu:settings=3",
       "comment": "Disable presolve and propagation for SCIP for duals."
     }
   },
diff --git a/test/end2end/scripts/python/Solver.py b/test/end2end/scripts/python/Solver.py
index 497227050..8d2c92a10 100644
--- a/test/end2end/scripts/python/Solver.py
+++ b/test/end2end/scripts/python/Solver.py
@@ -758,6 +758,23 @@ def __init__(self, exeName, timeout=None, nthreads=None, otherOptions=None):
         super().__init__(exeName, timeout, nthreads, otherOptions, stags)
 
 
+class GCGSolver(MPDirectSolver):
+    #def _setLPMethod(self, method : str):
+    #    m  = "simplex" if method == "SIMPLEX" else "ipm"
+    #    return f"alg:method {m}"
+
+    def _getAMPLOptionsName(self):
+        return "gcg"
+
+    def _setNThreads(self, threads):
+        return ""
+
+    def __init__(self, exeName, timeout=None, nthreads=None, otherOptions=None):
+        stags = {ModelTags.continuous, ModelTags.integer, ModelTags.binary,
+                 ModelTags.return_mipgap}
+        super().__init__(exeName, timeout, nthreads, otherOptions, stags)
+
+
 class HighsSolver(MPDirectSolver):
     def _setLPMethod(self, method : str):
         m  = "simplex" if method == "SIMPLEX" else "ipm"
diff --git a/test/end2end/scripts/python/SolverCollection.py b/test/end2end/scripts/python/SolverCollection.py
index 7302b12e6..dd9b9e1fe 100644
--- a/test/end2end/scripts/python/SolverCollection.py
+++ b/test/end2end/scripts/python/SolverCollection.py
@@ -41,6 +41,8 @@ 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.GCGSolver(path.join(binPath, "gcg")))
+    solvers.addSolver(Solver.GCGSolver(path.join(binPath, "gcgmp")))
     solvers.addSolver(Solver.SCIPSolver(path.join(binPath, "scip")))
     solvers.addSolver(Solver.SCIPSolver(path.join(binPath, "scipmp")))
     solvers.addSolver(Solver.CPLEXODHSolver(path.join(binPath, "cplexodh")))