/* 3party.c */

#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 */