mt_SL.c

Go to the documentation of this file.

/* MTgomory "multi tableau gomory cut" provides an implementation for gomory
 * cuts derived from multiple tableau rows in the spirit of the work of
 * Andersen et al (IPCO 2007), Cornuejols (es presented in George Nemhauser
 * Birthday Conference in Atlanta 2007) and Gomory (presented in the same
 * conference).
 *
 * Copyright (C) 2007 Daniel Espinoza.
 * 
 * This library is free software; you can redistribute it and/or modify it
 * under the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation; either version 2.1 of the License, or (at your
 * option) any later version.
 *
 * This library is distributed in the hope that it will be useful, but 
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public 
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this library; if not, write to the Free Software Foundation,
 * Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA 
 * */
/* ========================================================================= */
#include "mt_SL.h"
#include "mt_cplex_cbk.h"
/** @file 
 * @ingroup MTgomory */
/** @addtogroup MTgomory */
/** @{ */
/* ========================================================================= */
/** @brief compute in work the value ineq * r_i as described by the given
 * inequality and tableau, we assume that work is zero, moreover, given the
 * variable status (interger/binary) and an operation mode, try to do nothing
 * (mode 0), try to get ineq*r[i] >0 (mode 1) or try to get ineq*r[i] < 0 (mode
 * 2), note that mode != 0 will modify the tableau
 * @note We assume that the coefficients for integer variables in the tableau
 * are between ]-1:1[
 * @note when computing a new inequality, mode can be non-zero only at the
 * beggining, otherwise, the resulting inequality might be invalid.
 * @return zero on success, non-zero otherwise */
int MTsl1(
    const double*const ineq,
    MTrowmatrix_t*const tb,
    double*const work,
    const char*const vtype,
    const int mode)
{
  const int*const rowbeg = tb->rowbeg;
  const int*const rowind = tb->rowind;
  double*const rowval = tb->rowval;
  register int i, j;
  double dtmp;
  int col, isint, rval = 0;
  /* first tableau row */
  for( j = 0 ; j < tb->nrows ; j++)
  {
    dtmp = ineq[j];
    if(fabs(dtmp)<MT_ZERO_TOL) continue;
    for( i = rowbeg[j] ; i < rowbeg[j+1] ; i++)
    {
      col = rowind[i];
      work[col] += dtmp*rowval[i];
      isint = ((vtype[col]==CPX_BINARY) || (vtype[col]==CPX_INTEGER)) && (j==0);
      if(isint)
      {
        TESTGL( MT_DEBUG_LVL, (rval=(fabs(rowval[i])>1)), CLEANUP, "Imposible,"
                " coefficients in the tableau for integer variables should"
                " have |val = %lf| < 1", rowval[i]);
        if(mode > 0 && work[col] < 0)
        {
          work[col] -= dtmp*rowval[i];
          rowval[i] += (dtmp<0 ? -1.0 : 1.0);
          work[col] += dtmp*rowval[i];
        }
        else if( mode < 0 && work[col] > 0)
        {
          work[col] -= dtmp*rowval[i];
          rowval[i] += (dtmp<0 ? 1.0 : -1.0);
          work[col] += dtmp*rowval[i];
        }
      }/* end changing tableau according to signs */
    }/* end loop through non-zeros */
  }/* end loop through tableau rows */
  CLEANUP:
  return rval;
}
/* ========================================================================= */
/** @brief, set cutval[i] = max(cutval[i],work[i] > 0 ? work[i]/ineq[3]: 0);
 * and reset work to zero */
int MTsl2(
    const double*const ineq,
    double*const cutval,
    double*const work,
    const int * const active,
    const int nactive)
{
  register int i;
  int ccol, count=0;
  double dtmp;
  for( i = nactive ; i-- ; )
  {
    ccol = active[i];
    if(work[ccol]>0)
    {
      dtmp = work[ccol]/ineq[3];
      if(dtmp > cutval[ccol])
      {
        cutval[ccol] = dtmp;
        count++;
      }
    }
    work[ccol] = 0;
  }
  MESSAGE(MT_VERB_LVL,"Updated %d cut coeff", count);
  return 0;
}
/* ========================================================================= */
/** @brief given an inequality and f compute ineq[dim+1] and (possible change
 * signs) to ensure that ineq[dim+1] > 0; */
int MTsl3(
    const double*const f,
    double*const ineq,
    const int dim)
{
  register int i;
  int rval = 0;
  ineq[dim+1] = ineq[dim];
  for( i = dim ; i-- ; ) ineq[dim+1] -= ineq[i]*f[i];
  if(ineq[dim+1]<0)
  {
    MESSAGE(MT_VERB_LVL,"Change sign");
    for( i = dim+2; i-- ; ) ineq[i] = - ineq[i];
  }
  TESTG((rval=(fabs(ineq[dim+1])<MT_ZERO_TOL)), CLEANUP, 
        "f satisfy ax <= b at equality, a=[%lf,%lf] b = %lf f=[%lf,%lf]",
        ineq[0], ineq[1], ineq[dim], f[0], f[1]);
  CLEANUP:
  return rval;
}
/* ========================================================================= */
/** @brief compute the line that pases through the two given points */
int MTsl4(
  double*const ineq,
  const double x0,
  const double _y0,
  const double x1,
  const double _y1)
{
  const double Dx = (x0 - x1);
  const double Dy = (_y0 - _y1);
  ineq[0] = Dy;
  ineq[1] = -Dx;
  ineq[2] = x0*Dy - _y0*Dx;
  MESSAGE(MT_VERB_LVL,"x*%lg + y*%lg = %lg", ineq[0], ineq[1], ineq[2]);
  return 0;
}
/* ========================================================================= */
/** @brief compute active set from the given compresed 
 * vector set */
int MTsl5(
  int*const nactive,
  int*const active,
  const int nvars,
  const int nz,
  const double*const compval,
  const int*const compind)
{
  register int i;
  int cnt;
  memset(active,0,sizeof(int)*((size_t)nvars));
  for( i = nz ; i-- ; ) if(fabs(compval[i])>MT_ZERO_TOL) active[compind[i]] = 1;
  for(cnt = 0, i=0 ; i < nvars ; i++ ) if(active[i]) active[cnt++] = i;
  *nactive = cnt;
  return 0;
}
/* ========================================================================= */
int MTslCut(
    MTgomory_ccbk_t*const data,
    CPXCENVptr env,
    CPXLPptr lp)
{
  /* constant data */
  const MTsol_t*const sol = &(data->sol);
  const double*const x = sol->x;
  const double*const slack = sol->slack;
  const MTlp_t*const MTlp = &(data->lp);
  const int nrows = MTlp->nrows;
  const int ncols = MTlp->ncols;
  const int nvars = ncols + nrows;
  const char*const rtype = MTlp->rtype;
  const char*const ctype = MTlp->ctype;
  const int*const cstat = MTlp->cstat;
  MTcutHeap_t*const cuth = &(data->cuth);
  MTrowmatrix_t*const tb = &(data->tb);
  double*const f = data->f;
  double*const work = data->work;
  double*const cutval = data->extcut;
  /* local data */
  register int i, j;
  MTcut_t*cut = 0;
  int nfrac = 0;/* fractional integer variables in current solution*/
  int nint = 0;/* integer basic variables */
  int*frac_ind = EGsMalloc(int,ncols);/* fractional index set */
  int*int_ind = 0;  /* basic integral index set */
  double*binvrow = EGsMalloc(double, nrows);
  int rval = 0, icol, fcol, ccol, brow, *tabind = 0, status, tnz, fit, iit, onz;
  double*tableau = 0, ratio, minratio = DBL_MAX, dtmp, ineq[3][4], mval[2], util[2];
  int*active = EGsMalloc(int,nvars);
  char*vtype = EGsMalloc(char,nvars);
  int nactive = 0;
  /* build vtype */
  for( i = ncols ; i-- ; ) vtype[i] = ctype[i];
  vtype += ncols;
  for( i = nrows ; i-- ; ) vtype[i] = rtype[i];
  vtype -= ncols;
  /* identify integer fractional and integer variables */
  rval = MTcompute_fractional_vars( MTlp, x, frac_ind, &nfrac, 
                                    MT_CCBK_MIN_FRAC);
  CHECKRVALG(rval,CLEANUP);
  int_ind = frac_ind + nfrac;
  rval = MTcompute_integer_vars(MTlp, x, int_ind, &nint, 
                                MT_CCBK_MIN_FRAC);
  CHECKRVALG(rval,CLEANUP);
  /* allocate memory */
  MTrowmatrix_rsz(tb,nrows,0);
  tb->rowbeg[0] = 0;
  /* we double loop on fractional variables and in integer valued basic
   * variables tableau, and generate the two dimmensional lifted gomory cut */
  for( iit = nint; iit-- ; )
  {
    tb->nz = tb->nrows = 0;
    tb->ncols = nvars;
    icol = int_ind[iit];
    if(cstat[icol] != CPX_BASIC) continue;
    /* get asociated basic row */
    rval = CPXgetijrow(env, lp, -1, icol, &brow);
    MTcplexCHECKRVALG(env,rval,CLEANUP);
    rval = CPXbinvrow(env, lp, brow, binvrow);
    MTcplexCHECKRVALG(env,rval,CLEANUP);
    MTrowmatrix_rsz(tb,nrows,tb->ncols);
    tableau = tb->rowval;
    tabind = tb->rowind;
    rval = MTcplex_binv_to_tableau(MTlp, binvrow, tableau);
    CHECKRVALG(rval,CLEANUP);
    if(MT_VERB_LVL+10 < DEBUG) 
      MTcplex_display_tableau(stderr, MTlp, tableau);
    /* check that the asociated tableau is right */
    rval = MTcheckTabrow(MTlp,tableau,&status,&ratio,&tnz);
    CHECKRVALG(rval,CLEANUP);
    /* discard bad tableau rows */
    if(status) continue;
    if(minratio > ratio) minratio = ratio;
    /* compute f */
    MTcompute_f(MTlp, tableau, x, slack, f);
    /* now we compress the tableau */
    rval = MTcompressTabRow(MTlp, sol, tableau, tabind, f, &tnz);
    CHECKRVALG(rval,CLEANUP);
    MESSAGE(MT_VERB_LVL+5,"Tableau %d has %d non-zeros",icol,tnz);
    /* since these are integer tableau rows, they may end up being empty, in
     * such cases, we discard them */
    if(!tnz)
    {
      MESSAGE(MT_VERB_LVL,"Discarded empty integer tableau row %d",icol);
      continue;
    }
    /* up to now we have (but for integer factors) ax = f, and we want 
     * Z = f + Ax, moreover f[0] should be (about) zero */
    f[0] = round(f[0]) - f[0];
    if(fabs(f[0])>MT_CCBK_MIN_FRAC)
    {
      MTccbk_fail_tableau_frac++;
      MESSAGE(MT_VERB_LVL, "Discarding tableau %d, because f[%d] = %.4lf", 
          icol, 0, f[0]);
      continue;
    }
    tb->nrows = 1;
    onz = tb->nz = tnz;
    tb->rowbeg[1] = onz;
    /* now we loop through all fractional integer variables */
    for( fit = nfrac ; fit-- ; )
    {
      fcol = frac_ind[fit];
      /* check that the variable is basic */
      TESTG((rval=(cstat[fcol] != CPX_BASIC)), CLEANUP, "Variable[%d] = %lf,"
          " status %d isnt basic (%d)", fcol, x[fcol], cstat[fcol], CPX_BASIC); 
      /* get asociated basic row */
      rval = CPXgetijrow(env, lp, -1, fcol, &brow);
      MTcplexCHECKRVALG(env,rval,CLEANUP);
      rval = CPXbinvrow(env, lp, brow, binvrow);
      MTcplexCHECKRVALG(env,rval,CLEANUP);
      MTrowmatrix_rsz(tb,nrows,tb->nz+tb->ncols);
      tableau = tb->rowval+onz;
      tabind = tb->rowind+onz;
      rval = MTcplex_binv_to_tableau(MTlp, binvrow, tableau);
      CHECKRVALG(rval,CLEANUP);
      if(MT_VERB_LVL+10 < DEBUG) 
        MTcplex_display_tableau(stderr, MTlp, tableau);
      /* check that the asociated tableau is right */
      rval = MTcheckTabrow(MTlp,tableau,&status,&ratio,&tnz);
      CHECKRVALG(rval,CLEANUP);
      /* discard bad tableau rows */
      if(status) continue;
      if(minratio > ratio) minratio = ratio;
      /* store computed tableau */
      /* now we know that the tableau is good, and its length, so we save it */
      MTcompute_f(MTlp, tableau, x, slack, f+1);
      /* now we compress the tableau */
      rval = MTcompressTabRow(MTlp, sol, tableau, tabind, f+1, &tnz);
      CHECKRVALG(rval,CLEANUP);
      MESSAGE(MT_VERB_LVL+5,"Tableau %d has %d non-zeros",fcol,tnz);
      /* up to now we have (but for integer factors) ax = f, and we want 
       * Z = f + Ax */
      f[1] = -f[1] + round(f[1]);
      dtmp = fabs(f[1]);
      if(dtmp<MT_CCBK_MIN_FRAC)
      {
        MTccbk_fail_tableau_frac++;
        MESSAGE(MT_VERB_LVL, "Discarding tableau %d, because f[%d] = %.4lf", 
            fcol, 1, f[1]);
        continue;
      }
      /* f[1] should be between ]0,1[ */
      if(f[1]<0) f[1] += 1;
      /* now we passed all test for the tableau, so we save it */
      tb->nrows = 2;
      tb->nz = onz + tnz;
      tb->rowbeg[2] = tb->nz;
      #warning we could do a better job on integer variables in the second tableau row
      /* now we have the two row relaxation stored in tb, we now compute the 
       * cut */
      /* first see what active vars we have and reset work[active] == 0 */
      MESSAGE( MT_VERB_LVL, "joint tableau has %d rows and %d nonz", tb->nrows,
              tb->nz);
      rval = MTsl5(&nactive, active, nvars, tb->nz, tb->rowval, tb->rowind);
      CHECKRVALG(rval,CLEANUP);
      MESSAGE( MT_VERB_LVL, "Active cols %d %.3lf%%", nactive, 
              100*((double)nactive)/ncols);
      /* if nactive >= tb->nz, no intersection in tableau is possible */
      if(nactive >= tb->nz)
      {
        MESSAGE(MT_VERB_LVL,"Discarded, no intersection");
        continue;
      }
      #if MT_DEBUG_LVL +5 <= DEBUG
      FILE* fout = EGsfopen("cur_tableau.txt","w+");
      rval = MTwriteRTableau(fout, tb->nrows, tb->ncols, tb->rowval, tb->rowind, tb->rowbeg, f);
      CHECKRVALG(rval,CLEANUP);
      fclose(fout);
      #endif
      for( j = -1 ; j<2 ; j+=2 )
      {
        /* reset work and cutval */
        memset(work,0,sizeof(double)*((size_t)nvars));
        memset(cutval,0,sizeof(double)*((size_t)nvars));
        /* case j, set ineq[0] to j*x[0] <= 1 */
        /* set ineq[0] */
        MTsl4(ineq[0], (double)j, 1.0, (double)j, 0.0);
        /* b_i - a_i*f */
        rval = MTsl3(f,ineq[0],2);
        CHECKRVALG(rval,CLEANUP);
        /* compute in work[i] the value ineq[0] * r[i] */
        rval = MTsl1(ineq[0], tb, work, vtype, 1);
        CHECKRVALG(rval,CLEANUP);
        /* re-compute active */
        rval = MTsl5(&nactive, active, nvars, tb->nz, tb->rowval, tb->rowind);
        CHECKRVALG(rval,CLEANUP);
        MESSAGE(MT_VERB_LVL,"Active cols %d %.3lf%%", nactive, 
                100*((double)nactive)/ncols);
        /* compute current cutval and reset work */
        rval = MTsl2(ineq[0], cutval, work, active, nactive);
        CHECKRVALG(rval,CLEANUP);
        if(MT_VERB_LVL+5<=DEBUG) MTcplex_display_tableau(stderr,MTlp,cutval);
        /* compute mval */
        mval[0] = 1;
        mval[1] = 0;
        util[0] = 0;
        util[1] = 1;
        /* compute in work[i] the value util * r[i], and resero work to zero */
        rval = MTsl1(util, tb, work, vtype, 0);
        CHECKRVALG(rval,CLEANUP);
        for( i = nactive ; i-- ; )
        {
          ccol = active[i];
          dtmp = work[ccol]*cutval[ccol];
          if(dtmp > 0)
          {
            mval[0] = EGmax(mval[0],dtmp+f[1]);
          }
          else
          {
            mval[1] = EGmin(mval[1],dtmp+f[1]);
          }
          work[ccol] = 0;
        }
        /* if Delta mval is small, discard lifting */
        #if 1
        dtmp = mval[0] - mval[1];
        if(dtmp < 1.01)
        {
          MESSAGE(MT_VERB_LVL, "Discarding lifting Delta mval = %lg < 1.1", dtmp);
          continue;
        }
        #endif
        /* now we can compute ineq[1] as passing through (j,mval[0]) 
         * and (0,1), while ineq[2] as passing through (j,mval[1]) 
         * and (0,0) */
        MESSAGE(MT_VERB_LVL, "computing ineq for points (%lf,%lf) to (%lf,%lf)",
                (double)j, mval[0], 0.0, 1.0);
        MTsl4(ineq[1], (double)j, mval[0], 0.0, 1.0);
        /* b_i - a_i*f */
        rval = MTsl3(f,ineq[1],2);
        CHECKRVALG(rval,CLEANUP);
        /* compute in work[i] the value ineq[0] * r[i] */
        rval = MTsl1(ineq[1], tb, work, vtype, 0);
        CHECKRVALG(rval,CLEANUP);
        /* compute current cutval and reset work */
        rval = MTsl2(ineq[1], cutval, work, active, nactive);
        CHECKRVALG(rval,CLEANUP);
        if(MT_VERB_LVL+5<=DEBUG) MTcplex_display_tableau(stderr,MTlp,cutval);
        MESSAGE(MT_VERB_LVL, "computing ineq for points (%lf,%lf) to (%lf,%lf)",
                (double)j, mval[1], 0.0, 0.0);
        MTsl4(ineq[2], (double)j, mval[1], 0.0, 0.0);
        /* b_i - a_i*f */
        rval = MTsl3(f,ineq[2],2);
        CHECKRVALG(rval,CLEANUP);
        /* compute in work[i] the value ineq[0] * r[i] */
        rval = MTsl1(ineq[2], tb, work, vtype, 0);
        CHECKRVALG(rval,CLEANUP);
        /* compute current cutval and reset work */
        rval = MTsl2(ineq[2], cutval, work, active, nactive);
        CHECKRVALG(rval,CLEANUP);
        if(MT_VERB_LVL+5<=DEBUG) MTcplex_display_tableau(stderr,MTlp,cutval);
        /* now we must process the cut at hand */
        cut = MTcutHeap_get_free_cut(cuth);
        rval = MTraw_to_uncomplemented_cut( MTlp, cutval, nactive, active, 
                                            cut, &status);
        CHECKRVALG(rval,CLEANUP);
        TESTG((rval=(cut->nz != nactive)), CLEANUP, 
              "Imposible, cutnz = %d active %d", cut->nz, nactive);
        /* now we add the cut to the cutheap */
        if(status == 0)
        {
          rval = MTcutHeap_add_cut(MTlp,cuth);
          CHECKRVALG(rval,CLEANUP);
        }
      }/* end looping through both possible simple siftings */
    }/* end loop through set of fractional integer variables */
  }/* end loop through the set of integer vaariables */
  /* initialize values */
  CLEANUP:
  EGfree(frac_ind);
  EGfree(binvrow);
  EGfree(vtype);
  EGfree(active);
  return rval;
}
/* ========================================================================= */
int MTsl_ccbk(CPXCENVptr env,
              void*cbdata,
              int wherefrom,
              void*cbhandle,
              int*useraction_p)
{
  MTgomory_ccbk_t*const data=(MTgomory_ccbk_t*)cbhandle;
  MTcutHeap_t*const cuth = &(data->cuth);
  MTcut_t*cut=0;
  int rval = 0, action = 0;
  CPXLPptr lp = 0;
  int ncols,nrows; 
  /* set default action */ 
  MESSAGE(MT_VERB_LVL,"Entering");
  *useraction_p = CPX_CALLBACK_DEFAULT; 
  /* process callback info */
  rval = MTccbk_info_process(&(data->info), env, cbdata, wherefrom, &action);
  CHECKRVALG(rval,CLEANUP);
  if(action) goto CLEANUP;
  /* get the LP */
  rval = CPXgetcallbacknodelp (env, cbdata, wherefrom, &lp); 
  MTcplexCHECKRVALG(env,rval,CLEANUP);
  /* if we are verbose, write the problem to a file */
  if(MT_VERB_LVL +5 <= DEBUG)
  {
    rval = CPXwriteprob(env,lp,"node_lp.lp.gz","LP");
    MTcplexCHECKRVALG(env,rval,CLEANUP);
  }
  /* and get its description */
  rval = MTlp_load_problem(env, &(data->lp), lp, cbdata, wherefrom);
  CHECKRVALG(rval,CLEANUP);
  /* get numrows and numcols */
  ncols = CPXgetnumcols(env,lp);
  nrows = CPXgetnumrows(env,lp);
  MESSAGE(MT_VERB_LVL+3,"Problem has %d vars and %d constr",ncols,nrows);
  /* resize */
  rval = EGmax(nrows,ncols);
  MTgomory_ccbk_rsz(data,nrows,ncols+nrows,rval);
  /* now we get the current X and slack */
  rval = MTget_sol(&(data->lp), env, lp, cbdata, wherefrom, &(data->sol));
  CHECKRVALG(rval,CLEANUP);
  /* if we are debugging, display variable information */
  if(MT_VERB_LVL + 10 <= DEBUG)
    MTdisplay_lp(&(data->lp), &(data->sol), stderr);
  /* now we get the cuts */
  rval = MTslCut(data, env, lp);
  CHECKRVALG(rval,CLEANUP);
  /* now we add the cuts in the cutheap */
  data->info.added_cuts += MTcutHeap_get_ncuts(cuth);
  data->info.node_cuts += MTcutHeap_get_ncuts(cuth);
  /* now we are in shape to add the cut to CPLEX */
  while((cut=MTcutHeap_pop_cut(cuth)))
  {
    rval = MTuncomplemented_to_complemented_cut(&(data->lp), &(data->sol), cut);
    CHECKRVALG(rval,CLEANUP);
    WARNINGL(MT_VERB_LVL,!cut->nz,"empty cut!");
    if(!cut->nz) continue;
    /* test the cut extensivewlly if needed */
    if(MT_DEBUG_LVL + 10 <= DEBUG)
    {
      rval = MTcplex_test_cut(cut, &(data->lp));
      CHECKRVALG(rval,CLEANUP);
    }
    rval = CPXcutcallbackaddlocal( env, cbdata, wherefrom, cut->nz, cut->rhs, 
                              cut->sense , cut->cutind, cut->cutval);
    MTcplexCHECKRVALG(env,rval,CLEANUP);
    MESSAGE(MT_VERB_LVL, "Node %zd LP:(%d,%d,%d) Constr %p: %d nz, %.2lf rhs,"
          " ratio %.2lf, violation (>0) %lf", data->info.node_cnt, 
          data->lp.nrows, data->lp.ncols, data->lp.nz, (void*)cut, cut->nz,
          cut->rhs, cut->max_abs/cut->min_abs, cut->violation);
  }
  CLEANUP:
  MESSAGE(MT_VERB_LVL,"Done, status %d",rval);
  return rval;
}
/** @} */
/* end of mt_SL.c */
/* ========================================================================= */

---