mt_T2.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_T2.h"
/** @file 
 * @ingroup MTgomory */
/** @addtogroup MTgomory */
/** @{ */
/* ========================================================================= */
int MTt2Cut(const MTrowmatrix_t*const tb, 
            MTlp_t*const lp, 
            MTcutHeap_t*const cuth,
            const double*const f,
            double*const cutval,
            double*const work,
            const int cut_style)
{
  const int nrows = tb->nrows;
  const int ncols = tb->ncols;
  const double*const rowval = tb->rowval;
  const int*const rowind = tb->rowind;
  const int*const rowbeg = tb->rowbeg;
  int rval=0,pos,sense,redo,redo2, cur_tabrows = nrows;
  double* ineq_val = EGsMalloc(double,(ncols+1)*(nrows+1));
  double* ineq = 0;
  const int*order = 0;
  const int*sign = 0;
  int*active = EGsMalloc(int,ncols);
  int nactive = 0, status;
  MTcut_t*cut;
  register int i,j,k;
  EGgcIt_t gc;
  EGpermIt_t prit;
  double ratio = 0;
  #if MT_DEBUG_LVL +5 <= DEBUG
  FILE* fout = EGsfopen("cur_tableau.txt","w+");
  rval = MTwriteRTableau(fout, nrows, ncols, rowval, rowind, rowbeg, f);
  CHECKRVALG(rval,CLEANUP);
  fclose(fout);
  #endif
  /* check some basic assumptions */
  if(MT_DEBUG_LVL <= DEBUG)
  {
    for(i = nrows ; i-- ; )
    {
      TESTG((rval=(f[i]<=(-0x1p-1+MT_ZERO_TOL))), CLEANUP, 
            "%d-th fractional value is outside bounds %lf <= %lf", 
            i, f[i], (-0x1p-1+MT_ZERO_TOL));
      TESTG((rval=(f[i]>=(0x1p-1-MT_ZERO_TOL))), CLEANUP, 
            "%d-th fractional value is outside bounds %lf >= %lf", 
            i, f[i], (0x1p-1-MT_ZERO_TOL));
    }
  }
  /* allocate memory */
  REDO:
  EGgcItInit(&gc,cur_tabrows);
  EGpermItInit(&prit,cur_tabrows);
  order = EGpermItGetTuple(&prit);
  sign = EGgcItGetTuple(&gc);
  /* compute number of active columns in the tableau */
  memset(active,0,sizeof(int)*((size_t)ncols));
  for( i = rowbeg[cur_tabrows] ; i-- ;)
    active[rowind[i]] = 1;
  for(nactive=0, i = 0 ; i < ncols ; i++)
    if(active[i]) active[nactive++] = i;
  MESSAGE(MT_VERB_LVL,"Active cols %.3lf%%",100*((double)nactive)/ncols);
  /* man iterations */
  do
  {
    /* re-initialize gc iterator */
    EGgcItReset(&gc);
    if(MT_DEBUG_LVL + 10 <= DEBUG)
    {
      fprintf(stderr,"\rUsing Order  ");
      for( i = 0 ; i<cur_tabrows ; i++ ) fprintf(stderr,"%d ",order[i]);
      fprintf(stderr,"\n");
    }
    /* now we have to initialize the value of each inequality for the current
     * permutation, position k*(ncols+1)+i , contains 
     * const[perm]_k*\sigma*col_i for i=0,ncols-1, and position 
     * k*(ncols+1)+ncols contains rhs_k - const[perm]_k*\sigma*f */
    memset(ineq_val,0,sizeof(double)*((size_t)((ncols+1)*(cur_tabrows+1))));
    /* initialize positions k*(ncols+1)+ncols */
    for( k = cur_tabrows ; k-- ; ) 
      ineq_val[k*(ncols+1)+ncols] = 0.5*((double)(k+1));
    ineq_val[cur_tabrows*(ncols+1)+ncols] = 0.5*((double)cur_tabrows);
    /* compute coefficients for each column and rhs in each constraint */
    for( j = cur_tabrows ; j-- ; )
    {
      sense = sign[j] ? -1:1;
      k = order[j];
      ineq_val[k*(ncols+1)+ncols] += sense*f[j];
      for( k++ ; k <= cur_tabrows ; k++)
        ineq_val[k*(ncols+1)+ncols] -= sense*f[j];
      for( i = rowbeg[j] ; i < rowbeg[j+1] ; i++)
      {
        k = order[j];
        ineq_val[k*(ncols+1)+rowind[i]] -= sense*rowval[i];
        for(k++ ; k <= cur_tabrows ; k++)
          ineq_val[k*(ncols+1)+rowind[i]] += sense*rowval[i];
      }
    }
    /* now we start computing inequalities for the current permutation in all
     * posible axis reorientations */
    do
    {
      if(MT_DEBUG_LVL + 10 <= DEBUG)
      {
        fprintf(stderr,"\tUsing orientation ");
        for( i = 0 ; i<cur_tabrows ; i++ ) fprintf(stderr,"%d ",sign[i]);
        fprintf(stderr,"\n");
      }
      /* we have all values, so we can compute the resulting cut */
      /* compute the current cut in cutval */
      memset(cutval,0,sizeof(double)*((size_t)ncols));
      for( k = cur_tabrows+1 ; k-- ; )
      {
        ineq = ineq_val+k*(ncols+1);
        TESTG((rval = (ineq[ncols]<0)), CLEANUP, 
              "Imposible rhs_k - const[perm]_k*\\sigma*f < 0!");
        for( i = nactive ; i-- ; ) 
        {
          pos = active[i];
          if(ineq[pos] < MT_ZERO_TOL) continue;
          ratio = ineq[pos]/ineq[ncols];
          if(cutval[pos] < ratio) cutval[pos] = ratio;
        }
      }
      /* at this stage we have the cut computed and ineq_val, so we just
       * process the cut at hand. */
      cut = MTcutHeap_get_free_cut(cuth);
      for( i = ncols ; i-- ; ) work[i] = cutval[i]; 
      rval = MTraw_to_uncomplemented_cut( lp, work, nactive, active, cut, &status);
      CHECKRVALG(rval,CLEANUP);
      if(status) goto NEXT;
      /* now we add the cut to the cutheap */
      rval = MTcutHeap_add_cut(lp,cuth);
      CHECKRVALG(rval,CLEANUP);
      /* now we prepare for the next reorientation */
      NEXT:
      /* test if there is a next reorientaation */
      redo = EGgcItNext(&gc);
      if(!redo) break;
      pos = EGgcItGetChange(&gc);
      sense = sign[pos] ? -1:1;
      /* update the values of ineq_val */
      k = order[pos];
      ineq_val[k*(ncols+1)+ncols] += 2*sense*f[pos];
      for( k++ ; k<= cur_tabrows ; k++)
        ineq_val[k*(ncols+1)+ncols] -= 2*sense*f[pos];
      for( i = rowbeg[pos]; i < rowbeg[pos+1] ; i++)
      {
        k = order[pos];
        ineq_val[k*(ncols+1)+rowind[i]] -= 2*sense*rowval[i];
        for(k++ ; k <= cur_tabrows ; k++)
          ineq_val[k*(ncols+1)+rowind[i]] += 2*sense*rowval[i];
      }
    } while(1);
    /* check if we have to do another permutation */
    redo2 = EGpermItNext(&prit);
    if(!redo2) break;
  } while(1);
  if(cut_style && cur_tabrows > 2)
  {
    EGgcItClear(&gc);
    EGpermItClear(&prit);
    cur_tabrows--;
    goto REDO;
  }
  /* and we are done! */
  CLEANUP:
  EGgcItClear(&gc);
  EGpermItClear(&prit);
  EGfree(ineq_val);
  EGfree(active);
  return rval;
}
/* ========================================================================= */
int MTt2_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;
  size_t ncpr=0; 
  /* 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);
  /* get the best MT_CCBK_MAX_VARS tableaus */
  rval = MTget_best_k_tableau_rows( &(data->lp), &(data->sol), 
                                    data->info.max_rows, env, lp, 
                                    &(data->tb), data->f);
  CHECKRVALG(rval,CLEANUP);
  if(!data->tb.nrows) goto CLEANUP;
  /* and now we should call the cut-generator routine */
  rval = MTt2Cut( &(data->tb), &(data->lp), cuth, data->f, 
                  data->extcut, data->work, data->info.cut_style);
  CHECKRVALG(rval,CLEANUP);
  /* now we add the cuts in the cutheap */
  ncpr = MTcutHeap_get_ncuts(cuth);
  data->info.added_cuts += ncpr;
  data->info.node_cuts += ncpr;
  /* 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);
    /* test the cut extensivewlly if needed */
      if(MT_DEBUG_LVL + 10 <= DEBUG)
    {
      rval = MTcplex_test_cut(cut, &(data->lp));
      if(rval)
      {
        fprintf(stderr,"Failing cut was:");
        MTcplex_display_compress_tableau( stderr, &(data->lp), cut->cutval, 
                                          cut->cutind, cut->nz);
        fprintf(stderr," >= %.3lf, max_abs %lf, min_abs %lf\n",cut->rhs, 
                cut->max_abs, cut->min_abs);
      }
      CHECKRVALG(rval,CLEANUP);
      fprintf(stderr,"Adding cut with raw score %.7lg",cut->score);
      fprintf(stderr," MTcut: ");
      MTcplex_display_compress_tableau(stderr, &(data->lp), cut->cutval, cut->cutind, cut->nz);
      fprintf(stderr," >= %.4lf\n", cut->rhs);
    }
    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: %d nz, %.2lf rhs,"
          " ratio %.2lf, violation (>0) %lf", data->info.node_cnt, 
          data->lp.nrows, data->lp.ncols, data->lp.nz, cut->nz, cut->rhs, 
          cut->max_abs/cut->min_abs, cut->violation);
  }
  CLEANUP:
  MESSAGE(MT_VERB_LVL,"Done, status %d, added %zd cuts",rval,ncpr);
  return rval;
}
/** @} */
/* end of mt_T2.c */
/* ========================================================================= */

---