GLPK/Mixing GLPK with other solver packages

From Wikibooks, open books for an open world
Jump to: navigation, search

Users, on occasion, prefer to generate their problems with GLPK, but then solve them with another solver. Or conversely, import problems into GLPK which have been generated by another solver package. This page concentrates on the options involving high-level programming.


Introduction[edit]

First, the term foreign solver package refers to products like CPLEX or lp_solve. The motivations for combining two packages are various and not discussed here.

The options for shuffling between GLPK and another solver package can be tackled on three levels:

  • shell scripting to invoke each solver in turn via its command-line, while exchanging data via the filesystem using a common problem format (like MPS)
  • high-level programming using API calls from both GLPK and the foreign package, but still exchanging data via the filesystem
  • high-level programming using API calls from both GLPK and the foreign package, but using in-memory data structures instead

More information on standard problem formats can be found on the interoperability page.

COIN-OR Open solver interface (Osi)[edit]

The Osi initiative is described elsewhere in this book.

Compiling the COIN-OR CBC solver with mathprog[edit]

Using versions 2.7 or later, COIN-OR CBC solver can compile in the GLPK mathprog language as an input modelling language, although at this time no output features, such as database output or console printf statements, are supported. The COIN-OR CBC solver has implemented multi-threaded operation for the MIP solution techniques and is currently able to solve much larger and more complex MIP problems than GLPK.

Version 2.7 of the COIN-OR CBC solver currently uses GLPK 4.45 although it is planned to use 4.47 for the 2.8 release.

To compile and run the CBC solver on linux or Max OSX platforms follow the following steps (stable release 2.7 using a fresh install of Ubuntu 12.10 64 bit as an example):

 ~$ sudo apt-get install build-essential
 ~$ sudo apt-get install autoconf
 ~$ sudo apt-get install subversion
 ~$ svn co https://projects.coin-or.org/svn/Cbc/stable/2.7 coin-Cbc
 ~$ cd coin-Cbc/ThirdParty/Glpk
 ~/coin-Cbc/ThirdParty/Glpk$ ./get.Glpk
 ~/coin-Cbc/ThirdParty/Glpk$ ./configure
 ~/coin-Cbc/ThirdParty/Glpk$ cd ../..
 ~/coin-Cbc$ ./configure CFLAGS='-m64 -O3' CC=gcc-4.7 --enable-gnu-packages --enable-cbc-parallel --enable-debug -C --prefix=/usr/local
 ~/coin-Cbc$ make -j 8 [if you want to use parallel compilation - just 'make' otherwise]
 ~/coin-Cbc$ sudo make install     

Now copy your mathprog model to the machine and run the model, writing all of the computed variables to a csv file:

 cbc modelname.mod% -randomCbcSeed 0 -randomSeed 0 -proximity on -threads 8 -printi csv -solve -solu modelname.csv

Note the '%' at the end of the model name indicating that it is a mathprog model.

Filesystem interchange[edit]

The simplest high-level method is to use the filesystem to transfer the problem between GLPK and the foreign solver, and then transfer the solution back

For example, GLPK could write out the problem in free MPS format format, call the external solver, and then load the solution back into a new (second) glpk problem object. In this case, the API routines glp_write_prob, glp_read_mip, and glp_write_mip should be useful. In the unlikely event that hard-drive I/O contributes to poor performance, a RAM-drive can be deployed.

In-memory method[edit]

If the foreign solver has a suitable API, then the entire exercise can take place in the program space.

The following example uses CPLEX. CPLEX is a leading commercial solver suite, purchased in 2009 by IBM. Here GLPK is used to build a GLPK problem, transform that problem into a CPLEX problem, call CPLEX to solve the problem, and finally recover the solution. The following public domain C code titled 3party.c does just this. Users should be cautioned that 3party.c uses internal (non-public) data structures from GLPK — although the code could be easily rewritten to use published API data structures. It is always bad practice to stray from the public interface.

3party.c
/* Filename : 3party.c
 * Date     : July 2011
 * Author   : Andrew Makhorin <mao@gnu.org>
 * Tested   : GLPK v4.45
 *
 * Waiver: To the extent possible under law, Andrew Makhorin
 * <mao@gnu.org> has waived all copyright and related or
 * neighboring rights to this program code.
 * http://creativecommons.org/publicdomain/zero/1.0/
 *
 * Caution: This code uses internal (non-public) data structures
 * and should be rewritten to make use of published APIs instead.
 */
 
#include "cplex.h"
#include "glpapi.h"
 
static void cplex_error(CPXENVptr env, const char *who, int status)
{     char buffer[511+1];
      if (CPXgeterrorstring(env, status, buffer) == NULL)
         xprintf("%s failed\nCPLEX Error %5d: Unknown error code.\n",
            who, status);
      else
         xprintf("%s failed\n%s", who, buffer);
      return;
}
 
int solve_mip(glp_prob *P, const glp_iocp *parm)
{     CPXENVptr env;
      CPXLPptr prob = NULL;
      CPXFILEptr logfile = NULL;
      GLPROW *row;
      GLPCOL *col;
      GLPAIJ *aij;
      int status, numrows, numcols, objsen, numnz, m, n;
      int *matbeg, *matcnt, *matind, loc, i, j, k;
      char *sense, *ctype;
      int *ref = NULL, ret, *indices;
      double *matval, *obj, *lb, *ub, *rhs, *rngval, *x, sum;
      /* initialize CPLEX environment */
      env = CPXopenCPLEX(&status);
      if (env == NULL)
      {  xprintf("CPXopenCPLEX failed; CPLEX Error %5d.\n", status);
         ret = GLP_EFAIL;
         goto done;
      }
      xprintf("CPLEX version is %s\n", CPXversion(env));
      /* set CPLEX log file */
      logfile = CPXfopen("CONOUT$", "w");
      xassert(logfile != NULL);
      status = CPXsetlogfile(env, logfile);
      if (status)
      {  cplex_error(env, "CPXsetlogfile", status);
         ret = GLP_EFAIL;
         goto done;
      }
      /* create CPLEX problem object */
      prob = CPXcreateprob(env, &status, "MIP");
      if (prob == NULL)
      {  cplex_error(env, "CPXcreateprob", status);
         ret = GLP_EFAIL;
         goto done;
      }
      /* determine original number of rows and columns */
      m = P->m, n = P->n;
      /* build the row reference array;
         ref[i], i = 1,...,m, is the number which i-th row would have
         in the problem object after removing all free rows;
         if i-th row is free, ref[i] is set to 0;
         also count the number of rows and constraint coefficients in
         CPLEX problem object */
      ref = xcalloc(1+m, sizeof(int));
      numrows = 0, numnz = 0;
      for (i = 1; i <= m; i++)
      {  row = P->row[i];
         if (row->type == GLP_FR)
            ref[i] = 0;
         else
         {  ref[i] = ++numrows;
            for (aij = row->ptr; aij != NULL; aij = aij->r_next)
               numnz++;
         }
      }
      /* the set of columns includes one additional column fixed at 1
         to account the constant term of the objective function */
      numcols = n+1;
      /* allocate working arrays */
      obj = xcalloc(numcols, sizeof(double));
      sense = xcalloc(numrows, sizeof(char));
      rhs = xcalloc(numrows, sizeof(double));
      rngval = xcalloc(numrows, sizeof(double));
      matbeg = xcalloc(numcols, sizeof(int));
      matcnt = xcalloc(numcols, sizeof(int));
      matind = xcalloc(numnz, sizeof(int));
      matval = xcalloc(numnz, sizeof(double));
      lb = xcalloc(numcols, sizeof(double));
      ub = xcalloc(numcols, sizeof(double));
      /* set objective sense */
      if (P->dir == GLP_MIN)
         objsen = CPX_MIN;
      else if (P->dir == GLP_MAX)
         objsen = CPX_MAX;
      else
         xassert(P != P);
      /* set row attributes */
      for (k = 1; k <= m; k++)
      {  row = P->row[k];
         i = ref[k];
         if (i == 0) continue;
         if (row->type == GLP_LO)
         {  sense[i-1] = 'G';
            rhs[i-1] = row->lb;
            rngval[i-1] = 0.0;
         }
         else if (row->type == GLP_UP)
         {  sense[i-1] = 'L';
            rhs[i-1] = row->ub;
            rngval[i-1] = 0.0;
         }
         else if (row->type == GLP_DB)
         {  sense[i-1] = 'R';
            rhs[i-1] = row->lb;
            rngval[i-1] = row->ub - row->lb;
         }
         else if (row->type == GLP_FX)
         {  sense[i-1] = 'E';
            rhs[i-1] = row->lb;
            rngval[i-1] = 0.0;
         }
         else
            xassert(row != row);
      }
      /* set column attributes */
      for (j = 1; j <= n; j++)
      {  col = P->col[j];
         obj[j-1] = col->coef;
         if (col->type == GLP_FR)
         {  lb[j-1] = -CPX_INFBOUND;
            ub[j-1] = +CPX_INFBOUND;
         }
         else if (col->type == GLP_LO)
         {  lb[j-1] = col->lb;
            ub[j-1] = +CPX_INFBOUND;
         }
         else if (col->type == GLP_UP)
         {  lb[j-1] = -CPX_INFBOUND;
            ub[j-1] = col->ub;
         }
         else if (col->type == GLP_DB)
         {  lb[j-1] = col->lb;
            ub[j-1] = col->ub;
         }
         else if (col->type == GLP_FX)
            lb[j-1] = ub[j-1] = col->lb;
         else
            xassert(col != col);
      }
      obj[numcols-1] = P->c0;
      lb[numcols-1] = ub[numcols-1] = 1.0;
      /* build the constraint matrix in column-wise format */
      loc = 0;
      for (j = 1; j <= n; j++)
      {  col = P->col[j];
         matbeg[j-1] = loc;
         for (aij = col->ptr; aij != NULL; aij = aij->c_next)
         {  i = ref[aij->row->i];
            if (i != 0)
            {  matind[loc] = i-1;
               matval[loc] = aij->val;
               loc++;
            }
         }
         matcnt[j-1] = loc - matbeg[j-1];
      }
      matbeg[numcols-1] = loc;
      matcnt[numcols-1] = 0;
      xassert(loc == numnz);
      /* copy problem data to CPLEX problem object */
      status = CPXcopylp(env, prob, numcols, numrows, objsen, obj, rhs,
         sense, matbeg, matcnt, matind, matval, lb, ub, rngval);
      /* free working arrays */
      xfree(obj);
      xfree(sense);
      xfree(rhs);
      xfree(rngval);
      xfree(matbeg);
      xfree(matcnt);
      xfree(matind);
      xfree(matval);
      xfree(lb);
      xfree(ub);
      /* check if copying problem data is successful */
      if (status)
      {  cplex_error(env, "CPXcopylp", status);
         ret = GLP_EFAIL;
         goto done;
      }
      /* change problem type to MIP */
      status = CPXchgprobtype(env, prob, CPXPROB_MILP);
      if (status)
      {  cplex_error(env, "CPXchgprobtype", status);
         ret = GLP_EFAIL;
         goto done;
      }
      /* set column types */
      indices = xcalloc(numcols+1, sizeof(int));
      ctype = xcalloc(numcols+1, sizeof(char));
      for (j = 1; j <= n; j++)
      {  col = P->col[j];
         indices[j-1] = j-1;
         if (col->kind == GLP_CV)
            ctype[j-1] = 'C';
         else if (col->kind == GLP_IV)
            ctype[j-1] = 'I';
         else
            xassert(col != col);
      }
      indices[numcols-1] = numcols-1;
      ctype[numcols-1] = 'C';
      status = CPXchgctype(env, prob, numcols, indices, ctype);
      xfree(indices);
      xfree(ctype);
      if (status)
      {  cplex_error(env, "CPXchgctype", status);
         ret = GLP_EFAIL;
         goto done;
      }
      /* set MIP node log display level */
      status = CPXsetintparam(env, CPX_PARAM_MIPDISPLAY, 4);
      if (status)
      {  cplex_error(env, "CPXsetintparam(CPX_PARAM_MIPDISPLAY)",
            status);
         ret = GLP_EFAIL;
         goto done;
      }
      /* set solution time limit */
      if (parm->tm_lim < INT_MAX)
      {  status = CPXsetdblparam(env, CPX_PARAM_TILIM,
            (double)parm->tm_lim / 1000.0);
         if (status)
         {  cplex_error(env, "CPXsetdblparam(CPX_PARAM_TILIM)",
               status);
            ret = GLP_EFAIL;
            goto done;
         }
      }
      /* try to solve the problem */
      status = CPXmipopt(env, prob);
      if (status)
      {  cplex_error(env, "CPXmipopt", status);
         ret = GLP_EFAIL;
         goto done;
      }
      /* determine solution status */
      status = CPXgetstat(env, prob);
      switch (status)
      {  case CPXMIP_OPTIMAL:
         case CPXMIP_OPTIMAL_TOL:
            P->mip_stat = GLP_OPT;
            break;
         case CPXMIP_TIME_LIM_FEAS:
            P->mip_stat = GLP_FEAS;
            break;
         case CPXMIP_TIME_LIM_INFEAS:
            P->mip_stat = GLP_UNDEF;
            ret = GLP_ETMLIM;
            goto done;
         default:
            xprintf("CPXgetstat returned %d\n", status);
            ret = GLP_EFAIL;
            goto done;
      }
      /* obtain column values */
      x = xcalloc(numcols, sizeof(double));
      status = CPXgetmipx(env, prob, x, 0, numcols-1);
      if (status)
      {  cplex_error(env, "CPXgetmipx", status);
         xfree(x);
         P->mip_stat = GLP_UNDEF;
         ret = GLP_EFAIL;
         goto done;
      }
      for (j = 1; j <= n; j++)
      {  double t;
         col = P->col[j];
         if (col->kind == GLP_IV)
         {  t = floor(x[j-1] + 0.5);
            xassert(fabs(x[j-1] - t) <= 1e-3);
            x[j-1] = t;
         }
         col->mipx = x[j-1];
      }
      xfree(x);
      /* calculate objective value */
      sum = P->c0;
      for (j = 1; j <= n; j++)
      {  col = P->col[j];
         sum += col->coef * col->mipx;
      }
      P->mip_obj = sum;
      /* calculate row values */
      for (i = 1; i <= m; i++)
      {  row = P->row[i];
         sum = 0.0;
         for (aij = row->ptr; aij != NULL; aij = aij->r_next)
            sum += aij->val * aij->col->mipx;
         row->mipx = sum;
      }
      ret = 0;
done: if (ref != NULL)
         xfree(ref);
      if (prob != NULL)
         CPXfreeprob(env, &prob);
      if (logfile != NULL)
         fclose(logfile);
      if (env != NULL)
         CPXcloseCPLEX(&env);
      return ret;
}
 
/* eof */

CPLEX interoperation[edit]

CPLEX is a proprietary optimization sofware package from IBM. This mid-2011 thread discusses the converting of solutions between CPLEX and GLPK.

MiniSat CNF-SAT solver[edit]

From version 4.46, GLPK bundles the MiniSat CNF-SAT or conjunctive normal form satisfiability problem solver developed by Niklas Eén and Niklas Sörensson at the Chalmers University of Technology, Sweden. MiniSat's permissive MIT license allows it to be distributed as part of GLPK. The MiniSat wrapper call is named glp_minisat1.

Version 4.47 introduced the glp_intfeas1 API which offers a native implementation of the MiniSat code — please be aware that this function was marked "tentative" in this release. The GLPSOL option --minisat now redirects to this new solver.

A satisfiability (SAT) problem involves determining if a given boolean formula can evaluate TRUE through the suitable assignment of its variables. The conjunctive normal form (CNF) indicates restrictions on the way the boolean formula is expressed. Boolean formulas expressed in other ways are convertible to CNF.

Official documentation[edit]

The official GLPK documentation file doc/cnfsat.pdf contains a full reference for this feature. That material is therefore not repeated here. See obtaining GLPK.

Architecture[edit]

GLPK reads and writes CNF satisfiability problems in DIMACS format (see also GLPK format). Certain restrictions are placed on the problem specification so that it remains a valid CNF-SAT problem. These restrictions can be explicitly checked with the glp_check_cnfsat public call and are automatically confirmed when the GLPSOL --minisat option is deployed.

Within GLPK, CNF-SAT is considered a special case of MIP, where all the variables are binary and all the constraints are covering inequalities. That is, a particular CNF-SAT problem can be stored in a GLPK problem object of type glp_prob as if it were an MIP instance. This, in particular, means that a CNF-SAT problem can be solved with the glp_intopt solver, but that a call to glp_minisat1 or glp_infeas1, being problem-specific, should be much more efficient.

In these latter cases, GLPK translates its MIP into a suitable format, calls the appropriate solver which then attempts to identify a feasible solution, and finally translates the solution back to the original problem instance for subsequent analysis and reporting.