mt_cutpool.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_cutpool.h"
/** @file 
 * @ingroup MTgomory */
/** @addtogroup MTgomory */
/** @{ */
/* ========================================================================= */
void MTcutHeapInit(MTcutHeap_t*const cuth,size_t sz)
{
  register int i;
  MTcutHeapCn_t*chcn;
  memset(cuth,0,sizeof(MTcutHeap_t));
  i = ((int)(cuth->all_sz = sz+1));
  dbl_EGeHeapInit(&(cuth->heap));
  dbl_EGeHeapChangeD(&(cuth->heap),2);
  dbl_EGeHeapResize(&(cuth->heap),cuth->all_sz);
  cuth->all_cn = EGsMalloc(MTcutHeapCn_t,cuth->all_sz); 
  EGeListInit(&(cuth->free_cn));
  for( ; i-- ; ){
    MTcutHeapCnInit(cuth->all_cn+i);
    EGeListAddAfter(&(cuth->all_cn[i].lcn),&(cuth->free_cn));
    }
  chcn = EGcontainerOf(cuth->free_cn.next,MTcutHeapCn_t,lcn);
  cuth->next_free = &(chcn->cut);
  EGeListDel(cuth->free_cn.next);
  cuth->control.cut_score = MTcsDefault;
  return;
}
/* ========================================================================= */
void MTcutHeapClear(MTcutHeap_t*const cuth)
{
  register int i = (int)(cuth->all_sz);
  dbl_EGeHeapResize(&(cuth->heap),0);
  dbl_EGeHeapClear(&(cuth->heap));
  for(;i--;) MTcutHeapCnClear(cuth->all_cn+i);
  EGfree(cuth->all_cn);
  memset(cuth,0,sizeof(MTcutHeap_t));
  return;
}
/* ========================================================================= */
/** @brief given a #MTcutHeapCn_t compute the value asociated with the stored
 * cut and store it in #MTcutHeapCn_t::hcn, inside #dbl_EGeHeapCn_t::val, the
 * function used to evaluate depend on the configuration
 * for the cut heap.  The smallest the value, the worst the cut is.
 * @param chcn structure where we operate
 * @param cuth cutheap controling the cut heap connector 
 * @return zero on success, non-zero otherwise 
 * */
static int MTcutHeapCn_compute_val( MTcutHeap_t*const cuth,
                                    MTcutHeapCn_t*const chcn)
{
  register int i;
  double val = 0;
  int rval = 0;
  switch(cuth->control.cut_score)
  {
    case MTcsRAW:
      chcn->hcn.val = chcn->cut.score;
      break;
    case MTcsVIOL:
      chcn->hcn.val = chcn->cut.violation;
      break;
    case MTcsINVLINF:
      switch(chcn->cut.sense)
      {
        case 'L':
        case 'l':
          chcn->hcn.val = -chcn->cut.rhs/chcn->cut.max_abs;
          break;
        case 'G':
        case 'g':
          chcn->hcn.val = chcn->cut.rhs/chcn->cut.max_abs;
          break;
        default:
          FTESTG((rval=1),CLEANUP,"unexpected sense %c",chcn->cut.sense);
          break;
      }
      break;
    case MTcsINVL2:
      for( i = chcn->cut.nz ; i-- ; ) 
        val += chcn->cut.cutval[i]*chcn->cut.cutval[i];
      switch(chcn->cut.sense)
      {
        case 'L':
        case 'l':
          chcn->hcn.val = -chcn->cut.rhs*chcn->cut.rhs/val;
          break;
        case 'G':
        case 'g':
          chcn->hcn.val = chcn->cut.rhs*chcn->cut.rhs/val;
          break;
        default:
          FTESTG((rval=1),CLEANUP,"unexpected sense %c",chcn->cut.sense);
          break;
      }
      break;
    default:
      FTESTG((rval=1), CLEANUP, "Unknown function selection criteria %d", 
              cuth->control.cut_score);
      break;
  }
  CLEANUP:
  return rval;
}
/* ========================================================================= */
/** @brief permutation sort for integers */
void MTintPermSort( const int sz,
                    int*const perm,
                    const int*const elem)
{
  int i, j, temp;
  int t;
  if (sz <= 1)
    return;

  EGswap (perm[0], perm[(sz - 1) / 2], temp);
  i = 0;
  j = sz;
  t = elem[perm[0]];
  for (;;)
  {
    do
      i++;
    while (i < sz && (elem[perm[i]] < t));
    do
      j--;
    while (t < elem[perm[j]]);
    if (j < i)
      break;
    EGswap (perm[i], perm[j], temp);
  }
  EGswap (perm[0], perm[j], temp);
  MTintPermSort (j, perm, elem);
  MTintPermSort (sz - i, perm + i, elem);
}
/* ========================================================================= */
/** @brief internal memory for dominance check */
static int* MTcut1_nz = 0;
static int* MTcut2_nz = 0;
void MTcutNzAlloc(void) __attribute__ ((constructor));
void MTcutNzAlloc(void)
{
  MTcut1_nz = int_EGlpNumAllocArray(100);
  MTcut2_nz = int_EGlpNumAllocArray(100);
}
void MTcutNzFree(void) __attribute__ ((destructor));
void MTcutNzFree(void)
{
  int_EGlpNumFreeArray(MTcut1_nz);
  int_EGlpNumFreeArray(MTcut2_nz);
}
/* ========================================================================= */
/** @brief helper to set dominance */ 
#define __MT_DOM do{\
  if(v2 + MT_ZERO_TOL < v1) \
  {\
    if(domset && dom != -1) return 0;\
    domset=1;\
    dom=-1;\
    equal = 0; \
  }\
  else if( v1 + MT_ZERO_TOL < v2)\
  {\
    if(domset && dom != 1) return 0;\
    domset=1;\
    dom=1;\
    equal = 0;\
  }\
  }while(0)
/* ========================================================================= */
int MTcut_check_dominance(const MTcut_t*const cut1, const MTcut_t*const cut2)
{
  /* we work with ax >= rhs inequalities */
  const size_t c1nz = ((size_t)(cut1->nz));
  const size_t c2nz = ((size_t)(cut2->nz));
  const int sense1 = (cut1->sense == 'l' || cut1->sense == 'L') ? -1 : 1;
  const int sense2 = (cut2->sense == 'l' || cut2->sense == 'L') ? -1 : 1;
  int dom=0,domset=0,equal=1,p1,p2,c1,c2;
  register int i,j;
  double v1,v2;
  /* if either cut is an equation, we return 0 */
  if( cut1->sense == 'e' || cut1->sense == 'E' || 
      cut2->sense == 'e' || cut2->sense == 'E') return 0;
  /* allocate memory if necesary */
  if(__EGlpNumArraySize(MTcut1_nz) < c1nz)
    MTcut1_nz = int_EGlpNumReallocArray( &MTcut1_nz, c1nz);
  if(__EGlpNumArraySize(MTcut2_nz) < c2nz)
    MTcut2_nz = int_EGlpNumReallocArray( &MTcut2_nz, c2nz);
  /* sort the indices in the cut */ 
  for( i = cut1->nz ; i-- ; ) MTcut1_nz[i] = i;
  for( i = cut2->nz ; i-- ; ) MTcut2_nz[i] = i;
  MTintPermSort(cut1->nz, MTcut1_nz, cut1->cutind);
  MTintPermSort(cut2->nz, MTcut2_nz, cut2->cutind);
  /* first see if by checking RHS we can set a dominance direction */
  v1 = -1*sense1*cut1->rhs;
  v2 = -1*sense2*cut2->rhs;
  __MT_DOM;
  /* now we check the contents of both cuts */
  i = j = 0;
  while( i < cut1->nz || j < cut2->nz)
  {
    /* set the values to compare and advance in i and j */
    if( i < cut1->nz)
    {
      p1 = MTcut1_nz[i];
      c1 = cut1->cutind[p1];
      v1 = sense1*cut1->cutval[p1];
    }
    else
    {
      p1 = INT_MAX;
      c1 = INT_MAX;
      v1 = 0.0;
    }
    if( j < cut2->nz)
    {
      p2 = MTcut2_nz[j];
      c2 = cut2->cutind[p2];
      v2 = sense2*cut2->cutval[p2];
    }
    else
    {
      c2 = INT_MAX;
      p2 = INT_MAX;
      v2 = 0.0;
    }
    if(c1 < c2)
    {
      v2 = 0.0;
      j--;
    }
    else if(c1 > c2)
    {
      v1 = 0.0;
      i--;
    }
    i++;
    j++;
    /* compare the values */
    __MT_DOM;
  }
  /* ending */
  if(equal) dom=-1;
  return dom;
}
/* ========================================================================= */
void MTcut_copy(MTcut_t*const dest,const MTcut_t*const src)
{
  register int i = src->nz;
  dest->rhs = src->rhs;
  dest->violation = src->violation;
  dest->min_abs = src->min_abs;
  dest->max_abs = src->max_abs;
  dest->nz = src->nz;
  dest->sense = src->sense;
  dest->form = src->form;
  if(dest->nzspace < dest->nz)
    MTcut_rsz(dest,dest->nz);
  while(i--)
  {
    dest->cutind[i] = src->cutind[i];
    dest->cutval[i] = src->cutval[i];
  }
}
/* ========================================================================= */
int MTcutHeap_add_cut(const MTlp_t*const lp, MTcutHeap_t*const cuth)
{
  int rval = 0;
  register int i = ((int)(cuth->heap.sz));
  MTcutHeapCn_t*chcn = EGcontainerOf(cuth->next_free,MTcutHeapCn_t,cut);
  dbl_EGeHeapCn_t*ccn;
  MTcut_t *cut1 = &(chcn->cut),*cut2;
  /* now we compute the score for the cut, and add it to the heap */
  rval = MTcutHeapCn_compute_val(cuth,chcn);
  CHECKRVALG(rval, CLEANUP);
  /* check for dominated cuts */
  if(cuth->control.dominance_check)
  {
    for( ;i--;)
    {
      cut2 = &((EGcontainerOf(cuth->heap.cn[i], MTcutHeapCn_t,hcn))->cut);
      rval = MTcut_check_dominance(cut1,cut2);
      switch(rval)
      {
        case -1:
          MESSAGE(MT_VERB_LVL, "incomming cut is dominated, "
                  "score 1 %lf scrore 2 %lf",chcn->hcn.val, 
                  cuth->heap.cn[i]->val);
          rval = 0;
          goto CLEANUP;
          break;
        case 0:
          break;
        case 1:
          MTcut_copy(cut2,cut1);
          dbl_EGeHeapChangeVal(&(cuth->heap),cuth->heap.cn[i],chcn->hcn.val);
          MESSAGE(MT_VERB_LVL, "incomming cut dominate, "
                  "score 1 %lf scrore 2 %lf",chcn->hcn.val, 
                  cuth->heap.cn[i]->val);
          rval = 0;
          goto CLEANUP;
          break;
      }
      rval = 0;
    }
  }
  /* otherwise, we add it */
  rval = dbl_EGeHeapAdd(&(cuth->heap),&(chcn->hcn));
  CHECKRVALG(rval, CLEANUP);
  MESSAGE(MT_VERB_LVL, "Adding cut %p with score %lf",(void*)(&(chcn->cut)), 
          chcn->hcn.val);
  /* if we have no more empty cuts, we erase the worst from the heap */
  if(EGeListIsEmpty(&(cuth->free_cn)))
  {
    ccn = dbl_EGeHeapGetMin(&(cuth->heap));
    TESTG((rval = !ccn),CLEANUP,"Unexpected NULL pointer");
    dbl_EGeHeapDel(&(cuth->heap),ccn);
    chcn = EGcontainerOf(ccn,MTcutHeapCn_t,hcn);
    MESSAGE(MT_VERB_LVL, "Removing cut %p with score %lf",(void*)(&(chcn->cut)),
            ccn->val);
    EGeListAddAfter(&(chcn->lcn),&(cuth->free_cn));
    if((MT_VERB_LVL+3 <= DEBUG) && (&(chcn->cut) != cut1))
    {
      fprintf(stderr,"Adding new MTcut (compress+complemented): ");
      MTcplex_display_compress_tableau(stderr,lp,cut1->cutval, cut1->cutind, cut1->nz);
      fprintf(stderr,">=1\n");
    }
  }
  else if(MT_VERB_LVL+3 <= DEBUG)
  {
    fprintf(stderr,"Adding new MTcut (compress+complemented): ");
    MTcplex_display_compress_tableau(stderr,lp,cut1->cutval, cut1->cutind, cut1->nz);
    fprintf(stderr,">=1\n");
  }
  chcn = EGcontainerOf(cuth->free_cn.next,MTcutHeapCn_t,lcn);
  cuth->next_free = &(chcn->cut);
  EGeListDel(cuth->free_cn.next);
  CLEANUP:
  return rval;
}
/* ========================================================================= */
MTcut_t* MTcutHeap_pop_cut(MTcutHeap_t*const cuth)
{
  MTcutHeapCn_t*chcn=0;
  MTcut_t* cut=0;
  dbl_EGeHeapCn_t*ccn=0;
  ccn = dbl_EGeHeapGetMin(&(cuth->heap));
  if(ccn)
  {
    chcn = EGcontainerOf(ccn,MTcutHeapCn_t,hcn);
    cut = &(chcn->cut);
    dbl_EGeHeapDel(&(cuth->heap),ccn);
    EGeListAddAfter(&(chcn->lcn),&(cuth->free_cn));
  }
  return cut;
}
/* ========================================================================= */
void MTcut_rsz(MTcut_t*const MTcut,
              const int MTnsz)
{
  if(MTcut->nzspace < MTnsz)
  {
    MTcut->nzspace = MTnsz;
    MTcut->cutind = EGrealloc(MTcut->cutind, sizeof(int)*((size_t)MTnsz));
    MTcut->cutval = EGrealloc(MTcut->cutval, sizeof(double)*((size_t)MTnsz));
  }
  return;
}
/* ========================================================================= */
void MTcut_clear(MTcut_t*const MTcut)
{
  EGfree(MTcut->cutval);
  EGfree(MTcut->cutind);
  memset(MTcut,0,sizeof(MTcut_t));
  return;
}

/* ========================================================================= */
/** @} */
/* end of mt_cutpool.c */
/* ========================================================================= */

---