mt_T1.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_T1.h"
/** @file 
 * @ingroup MTgomory */
/** @addtogroup MTgomory */
/** @{ */
/* ========================================================================= */
int MTt1Cut(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, cur_tabrow = nrows;
  int*active = EGsMalloc(int,ncols);
  int nactive = 0, status;
  MTcut_t*cut = 0;
  register int i,j;
  EGgcIt_t gc;
  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
  memset(&gc,0,sizeof(EGgcIt_t));
  /* 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 */
  /* compute number of active columns in the tableau */
  memset(active,0,sizeof(int)*((size_t)ncols));
  for( i = rowbeg[nrows] ; 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);
  /* initialize values */
  REDO:
  EGgcItInit(&gc,cur_tabrow);
  for( i = ncols ; i-- ; ) work[i] = 0.0;
  work[ncols] = 0.5*cur_tabrow;
  /* now we compute in works[ncols] = n/2 - f*\sigma  and compute in work[i] 
   * col_i*\sigma */
  for( j = cur_tabrow ; j-- ; )
  {
    sense = ((EGgcItGetTuple(&gc))[j]);
    if(sense){
      work[ncols] -= f[j];
      for( i = rowbeg[j] ; i < rowbeg[j+1] ; i++)
        work[rowind[i]] += rowval[i];
    }
    else
    {
      work[ncols] += f[j];
      for( i = rowbeg[j] ; i < rowbeg[j+1] ; i++)
        work[rowind[i]] -= rowval[i];
    }
  }
  /* now we start computing the inequalities */ 
  do
  {
    TESTG((rval=(work[ncols]<0)), CLEANUP, "Imposible!, n/2-sigma*f"
          " (%.6lf) < 0",work[ncols]);
    /* compute the cut values */
    for( i = ncols ; i-- ; ) cutval[i] = work[i] > MT_ZERO_TOL ? work[i]/work[ncols]:0.0;
    for( j = cur_tabrow ; j-- ; )
    {
      sense = ((EGgcItGetTuple(&gc))[j]);
      if(sense)
      {
        for( i = rowbeg[j] ; i < rowbeg[j+1] ; i++)
        {
          if((rowval[i] < -MT_ZERO_TOL) && 
            ((ratio = rowval[i]/(-0.5-f[j])) > cutval[rowind[i]]))
            cutval[rowind[i]] = ratio;
        }
      }
      else
      {
        for( i = rowbeg[j] ; i < rowbeg[j+1] ; i++)
        {
          if((rowval[i] > MT_ZERO_TOL) &&
            ((ratio = rowval[i]/(0.5-f[j])) > cutval[rowind[i]]))
            cutval[rowind[i]] = ratio;
        }
      }
    }
    /* now we must process the cut at hand */
    cut = MTcutHeap_get_free_cut(cuth);
    rval = MTraw_to_uncomplemented_cut( lp, cutval, 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:
    redo = EGgcItNext(&gc);
    if(!redo) break;
    /* now we should re-compute work[*] */
    pos = EGgcItGetChange(&gc);
    sense = ((EGgcItGetTuple(&gc))[pos]);
    if(sense)
    {
      work[ncols] -= 2*f[pos];
      for( i = rowbeg[pos] ; i < rowbeg[pos+1] ; i++)
        work[rowind[i]] += 2*rowval[i];
    }
    else
    {
      work[ncols] += 2*f[pos];
      for( i = rowbeg[pos] ; i < rowbeg[pos+1] ; i++)
        work[rowind[i]] -= 2*rowval[i];
    }
  }while(1);
  EGgcItClear(&gc);
  if(cut_style && cur_tabrow > 2)
  {
    cur_tabrow--;
    EGgcItClear(&gc);
    goto REDO;
  }
  CLEANUP:
  EGgcItClear(&gc);
  EGfree(active);
  return rval;
}
/* ========================================================================= */
int MTt1_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);
  /* 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 = MTt1Cut( &(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 */
  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);
    /* 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: %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",rval);
  return rval;
}
/** @} */
/* end of mt_T1.h */
/* ========================================================================= */

---