eg_greedykp.c

Go to the documentation of this file.

#include "eg_greedykp.h"
#define EG_PRIMALIZATION_ERROR 0.001

EGgreedyData_t* EGnewGreedyData(EGmemPool_t *mem)
{
   EGgreedyData_t *data;

   data = EGmemPoolSMalloc(mem, EGgreedyData_t, 1);

   data->bdG = 0;
   data->cycleG = 0;
   data->bdGnodeMap = 0;
   data->cycleEdgeMap = 0;
   data->ddom_list = 0;
   data->ee_dist = 0;
   data->ee_prec = 0;
   data->eo_prec = 0;

   return data;
}

void EGfreeGreedyData(void *vdata)
{

   unsigned int i;
   EGgreedyData_t *data = (EGgreedyData_t*)(vdata);

   EGmemPool_t *mem = data->cycleG->mem;
   if(data->randa) EGfree(data->randa);
   if(data->randb) EGfree(data->randb);

   if (data->pndummy)
      free(data->pndummy);
   if (data->pedummy)
      free(data->pedummy);

   if (data->bdnodes)
      free(data->bdnodes);
   if (data->pedges)
      free(data->pedges);
   if (data->pedgevals)
      free(data->pedgevals);

   if (data->bdGnodeMap)
      EGmemPoolSFree(data->bdGnodeMap, EGdGraphNode_t*, data->bdG->nodes->size, mem);
   if (data->cycleEdgeMap)
      EGmemPoolSFree(data->cycleEdgeMap, EGdGraphEdge_t*, data->cycleG->nedges, mem);

   for(i=0; i < data->bdG->nodes->size; i++)
   {
      if (data->ee_dist)
         EGmemPoolSFree(data->ee_dist[i], EGdijkstraCost_t, i+1, mem);
      #if EG_KP_HEURISTIC
      if (data->eo_dist)
         EGmemPoolSFree(data->eo_dist[i], EGdijkstraCost_t, i+1, mem);
      #endif
      if (data->ee_prec)
         EGmemPoolSFree(data->ee_prec[i], int, i+1, mem);
      if (data->eo_prec)
         EGmemPoolSFree(data->eo_prec[i], int, i+1, mem);
   }

   if (data->ee_dist)
      EGmemPoolSFree(data->ee_dist, EGdijkstraCost_t*, data->bdG->nodes->size, mem);
   #if EG_KP_HEURISTIC
   if (data->eo_dist)
      EGmemPoolSFree(data->eo_dist, EGdijkstraCost_t*, data->bdG->nodes->size, mem);
   #endif
   if (data->ee_prec)
      EGmemPoolSFree(data->ee_prec, int*, data->bdG->nodes->size, mem);
   if (data->eo_prec)
      EGmemPoolSFree(data->eo_prec, int*, data->bdG->nodes->size, mem);

   #if EG_KP_HEURISTIC
   EGmemPoolSFree(data->bdGtoCycleGmap, EGdGraphNode_t*, data->bdG->nodes->size, mem);
   #endif

   if (data->cycleG)
   {
      EGdGraphClearMP(data->cycleG, EGfreeCycleDataMP, EGfreeDijkstraNodeDataMP, 0, mem, mem, 0);
      EGfreeDGraph(data->cycleG);
   }

   EGmemPoolSFree(data, EGgreedyData_t, 1, mem);

   return;

}

int EGfillGreedyData (EGgreedyData_t * data,
                      EGdGraph_t * bdG,
                      EGlist_t * ddom_list,
                      EGlist_t** dembed,
                      graphP G,
                      int norig_nodes,
                      int norig_edges,
                      int nplanar_edges,
                      int *orig_edges,
                      int *planar_edges,
                      double *orig_weight,
                      double *planar_weight)
{

  int rval=0;
  register unsigned int i;
  EGmemPool_t *mem = bdG->mem;
  memset(data,0,sizeof(EGgreedyData_t));
  static size_t  dij_os[6] = { 
    [EG_DIJ_DIST] = offsetof (EGdijkstraNodeData_t, dist),
    [EG_DIJ_NDIST] = offsetof (EGdijkstraNodeData_t, ndist),
    [EG_DIJ_FATHER] = offsetof (EGdijkstraNodeData_t, father),
    [EG_DIJ_MARKER] = offsetof (EGdijkstraNodeData_t, marker),
    [EG_DIJ_HCONNECTOR] = offsetof (EGdijkstraNodeData_t, hc),
    [EG_DIJ_ELENGTH] = offsetof (EGcycleData_t, weight)
    };

  EGdijkstraCost_t max_cycle_edge = EGdijkstraToCost(1.0),
                   max_even_path = EGdijkstraToCost(1.0);

  data->norig_nodes = norig_nodes;
  data->norig_edges = norig_edges;
  data->nplanar_edges = nplanar_edges;
  data->randa = EGsMalloc(unsigned int,norig_edges);
  data->randb = EGsMalloc(unsigned int,norig_edges);
  data->pndummy = EGsMalloc(unsigned char, data->norig_nodes);
  data->pedummy = EGsMalloc(unsigned char, data->norig_edges);
  data->pedges = EGsMalloc(EGdGraphEdge_t*,data->norig_edges);
  data->bdnodes = EGsMalloc(EGdGraphNode_t*,bdG->nodes->size);
  data->pedgevals = EGsMalloc(double,data->norig_edges);

  data->nviolated = 0;
  data->min_violated =  10000.0;
  data->max_violated = -10000.0;
  
  for( i = norig_edges ; i-- ; )
  {
    data->randa[i] = random()>>12;
    data->randb[i] = random()>>12;
  }
  data->orig_edges = orig_edges;
  data->planar_edges = planar_edges;

  data->orig_weight = orig_weight;
  data->planar_weight = planar_weight;
 
  data->bdG = bdG;
  data->ddom_list = ddom_list;
  data->G = G;
  data->cycleG = EGnewCycleGraph (mem, ddom_list, bdG, max_cycle_edge);
  TEST(!data->cycleG, "error generating cycle graph");

#if DISPLAY_MODE
  fprintf (stdout, "KP greedy data: Cycle graph has %u nodes and %d arcs.\n",
                    data->cycleG->nodes->size, data->cycleG->nedges);
  fflush (stdout);
#endif

  data->dembed = dembed;
  data->ee_dist = EGmemPoolSMalloc (mem, EGdijkstraCost_t *, bdG->nodes->size);
  #if EG_KP_HEURISTIC
  data->eo_dist = EGmemPoolSMalloc (mem, EGdijkstraCost_t *, bdG->nodes->size);
  #endif
  data->ee_prec = EGmemPoolSMalloc (mem, int *, bdG->nodes->size);
  data->eo_prec = EGmemPoolSMalloc (mem, int *, bdG->nodes->size);

  for (i = 0; i < bdG->nodes->size; i++)
  {
    data->ee_dist[i] = EGmemPoolSMalloc (mem, EGdijkstraCost_t, i+1);
    #if EG_KP_HEURISTIC
    data->eo_dist[i] = EGmemPoolSMalloc (mem, EGdijkstraCost_t, i+1);
    #endif
    data->ee_prec[i] = EGmemPoolSMalloc (mem, int, i+1); 
    data->eo_prec[i] = EGmemPoolSMalloc (mem, int, i+1);
  }

  data->cycleEdgeMap = EGmemPoolSMalloc(mem, EGdGraphEdge_t*, data->cycleG->nedges);
  {
     EGlistNode_t *vit, *eit;
     EGdGraphNode_t *v;
     EGdGraphEdge_t *e;
     for( vit = data->cycleG->nodes->begin; vit; vit=vit->next )
     {
        v = (EGdGraphNode_t*)(vit->this);
        for(eit=v->out_edges->begin; eit; eit=eit->next)
        {
           e = (EGdGraphEdge_t*)(eit->this);
           data->cycleEdgeMap[e->id] = e;
        }
     }
  }

  /* Remember that the nodes in bdG may not be sorted according to their id.  *
   * Hence, it is necessary to use a map to establish which id is which.      */

  data->bdGnodeMap = EGmemPoolSMalloc(mem, EGdGraphNode_t*, bdG->nodes->size);
  {
     EGlistNode_t *it;

     for(i=0, it=bdG->nodes->begin; it; it=it->next, i++)
     {
        data->bdGnodeMap[i] = (EGdGraphNode_t*)(it->this);
     }
  }


#if EG_KP_HEURISTIC
  {
    EGdGraphNode_t *cycle_v, *bdg_v;
    EGlistNode_t *it;

    data->bdGtoCycleGmap = EGmemPoolSMalloc(mem, EGdGraphNode_t*, data->bdG->nodes->size);
    for( it = data->cycleG->nodes->begin; it; it=it->next->next )
    {
      cycle_v = (EGdGraphNode_t*)(it->this);
      bdg_v = data->bdGnodeMap[cycle_v->id / 2];
      data->bdGtoCycleGmap[ bdg_v->id ] = cycle_v;
    }
  }
#endif

  /* fill the data in the table */ 
#if EG_GREEDYKP_FILLTABLE
  {
     EGlistNode_t *u_it, *v_it;
     EGdGraphNode_t *s, *t;

     EGdijkstraCost_t dist;
     int prec;
     int s_id, t_id;

     EGheap_t *my_heap = EGnewHeap (mem, 2, data->cycleG->nodes->size);

     /* u_it will iterate through even nodes in cycleG */
     for(u_it=data->cycleG->nodes->begin; u_it; u_it=((u_it->next)->next))
     {

        TEST (!(u_it->next), "cycle graph badly defined");

        s = (EGdGraphNode_t*)(u_it->this);

        rval = EGpartialDijkstra (s,0,max_even_path,dij_os,my_heap,data->cycleG);
        TEST(rval, "error running dijkstra");

        /* v_it will iterate through nodes smaller than u_it in cycleG */
        for(v_it=data->cycleG->nodes->begin; v_it != u_it->next->next; v_it=v_it->next)
        {
           t = (EGdGraphNode_t*)(v_it->this);
          
           s_id = data->bdGnodeMap[s->id/2]->id;
           t_id = data->bdGnodeMap[t->id/2]->id;

           /* we should never request the even distance between a node and
            * itself! might as well make the software crash */

           if (v_it == u_it)
           {
              data->ee_dist[s_id][t_id] = 0;
              data->ee_prec[s_id][t_id] = -1;
              continue;
           }

           /* determine distance(s,t) and father(t) */ 
           if (EGdijkstraGetMarker(t, dij_os) == UINT_MAX)
           {
              dist = EG_DIJKSTRA_COST_MAX;
              prec = 0;
           }
           else
           {
              dist = EGdijkstraGetDist(t, dij_os);
              prec = EGdijkstraGetFather(t, dij_os)->id;
           }

           /* store distance(s,t) and father(t) in tables */
           if (t->id & 1)
           {
              if (t_id < s_id)
              {
                 #if EG_KP_HEURISTIC
                 data->eo_dist[s_id][t_id] = dist;
                 #endif
                 data->eo_prec[s_id][t_id] = prec;
              }
              else
              {
                 #if EG_KP_HEURISTIC
                 data->eo_dist[t_id][s_id] = dist;
                 #endif
                 data->eo_prec[t_id][s_id] = prec;
              }
           }
           else
           {
              if (t_id < s_id)
              {
                data->ee_dist[s_id][t_id] = dist;
                data->ee_prec[s_id][t_id] = prec;
              }
              else
              {
                data->ee_dist[t_id][s_id] = dist;
                data->ee_prec[t_id][s_id] = prec;
              }
           }
           
        }
     }

     EGfreeHeap (my_heap, mem);

  }
  /* done filling the data in the table */
#endif

  return rval;

}


EGdualCut_t* EGnewDualCut(EGmemPool_t *mem, unsigned int sz)
{

   EGdualCut_t *dc;

   dc = EGmemPoolSMalloc(mem, EGdualCut_t, 1);

   dc->sz = sz;
   dc->value = EGdijkstraToCost(0.0);

   if(sz) 
      dc->edges = EGmemPoolSMalloc(mem, EGdGraphEdge_t*, sz);
   else 
      dc->edges = 0;

   return dc;

}

void EGfreeDualCut(void *v, EGmemPool_t *mem)
{

   EGdualCut_t *dc = (EGdualCut_t*)(v);

   if (dc->sz)
      EGmemPoolSFree(dc->edges, EGdGraphEdge_t*, dc->sz, mem);
   EGmemPoolSFree(dc, EGdualCut_t, 1, mem);

   return;

}

EGdcutIter_t* EGnewDcutIter(EGgreedyData_t *gdata)
{

   EGdcutIter_t *zit;

   if (!gdata || !gdata->cycleG || !gdata->ddom_list)
      return 0;

   zit = EGmemPoolSMalloc(gdata->cycleG->mem, EGdcutIter_t, 1);

   zit->num_it = 0;
   zit->ddit = gdata->ddom_list->begin;
   zit->current_side = 0;
   zit->current_middle = 0;
   zit->gdata = gdata;

   return zit;

}

void EGfreeDcutIter(void *v)
{
 
   EGdcutIter_t *zit = (EGdcutIter_t*)(v);

   if (!v)
      return;

   EXIT(!zit->gdata || !zit->gdata->cycleG, "cycleG is null");

   EGmemPoolSFree(v, EGdcutIter_t, 1, zit->gdata->cycleG->mem);

   return;

}

int EGincrementDcutIter(EGdcutIter_t *zit)
{

   TEST(!zit || !zit->ddit, "error with zit");

   zit->num_it += 1;
   zit->ddit = zit->ddit->next;

   if (!zit->ddit)
   {
   
      if (!zit->current_side)
         zit->current_side = 1;
      else if (zit->current_middle < 2)
      {
         zit->current_middle += 1;
         zit->current_side = 0;
      }
      else 
         return 0;

      zit->ddit = zit->gdata->ddom_list->begin;

   }

   return 0;

}

EGdualCut_t* EGgetDcut(EGdcutIter_t *zit)
{

   unsigned int i, j, p1, p2;
   EGdualCut_t *dc;
   EGddomino_t *ddom;
   EGmemPool_t *mem = zit->gdata->cycleG->mem;

   size_t mos[6];
   mos[EG_DIJ_ELENGTH] = offsetof(EGmengerEdgeData_t, cost);

   EGdijkstraCost_t val=EGdijkstraToCost(0.0);

   EXIT(!zit || !zit->ddit, "bad zit");

   ddom = zit->ddit->this;

   switch(zit->current_middle)
   {
      case 0:
         p1=1; p2=2;
         break;
      case 1:
         p1=0; p2=2;
         break;
      case 2:
         p1=0; p2=1;
         break;
      default:
         return 0; 
         break;
   }

   dc = EGnewDualCut(mem, ddom->npath[p1]+ddom->npath[p2]);

   for(i=0, j=0; i < ddom->npath[p1]; i++, j++)
   {
      dc->edges[j] = (EGdGraphEdge_t*)(ddom->path[p1][i]);
      val = EGdijkstraCostAdd(val, EGdijkstraGetEdgeLength(dc->edges[j], mos));
   }
   for(i=ddom->npath[p2]; i; i--, j++)
   {
      dc->edges[j] = (EGdGraphEdge_t*)(ddom->path[p2][i-1]);
      val = EGdijkstraCostAdd(val, EGdijkstraGetEdgeLength(dc->edges[j], mos));
   }

   dc->value = val;

   return dc;

}

int EGgetOddPrec (unsigned int s, unsigned int t, EGgreedyData_t*data)
{

   if (t < s)
      return ( data->eo_prec[s][t] ); 
   else
      return ( data->eo_prec[t][s] ); 

}

int EGgetEvenPrec (unsigned int s, unsigned int t, EGgreedyData_t*data)
{
   #if DEBUG > 3 
   if (s >= data->bdG->nodes->size)
      MESSAGE(1, "out of range s %u / %u", s, data->bdG->nodes->size);
   if (t >= data->bdG->nodes->size)
      MESSAGE(1, "out of range t %u / %u", t, data->bdG->nodes->size);
   #endif
   
   if (t < s)
      return ( data->ee_prec[s][t] ); 
   else
      return ( data->ee_prec[t][s] ); 
}

EGdijkstraCost_t EGgetEvenDistance (unsigned s, unsigned t, EGgreedyData_t*data)
{

   #if DEBUG > 3
   if (s >= data->bdG->nodes->size)
      MESSAGE(1, "out of range s %u / %u", s, data->bdG->nodes->size);
   if (t >= data->bdG->nodes->size)
      MESSAGE(1, "out of range t %u / %u", t, data->bdG->nodes->size);
   #endif

   if(t < s)
      return ( data->ee_dist[s][t] ); 
   else
      return ( data->ee_dist[t][s] ); 

}

#if EG_KP_HEURISTIC
EGdijkstraCost_t EGgetOddDistance (unsigned s, unsigned t, EGgreedyData_t*data)
{

   if(t < s)
      return ( data->eo_dist[s][t] ); 
   else
      return ( data->eo_dist[t][s] ); 

}
#endif

int EGextractEEpath(unsigned int s, 
                    unsigned int t, 
                    EGgreedyData_t *data, 
                    EGlist_t *path)
{

   unsigned int curr_s = s, curr_t = t;
   int tot_even = 1, e_even = 1, rval=1;
   EGdGraphEdge_t *e;
   EGcycleData_t *e_data;

   EGdijkstraCost_t dist;
   dist = EGdijkstraToCost(0.0);

   EGlistClear(path, nullFree);

   int ncycle = 0;
   while ( curr_t != curr_s )
   {
    
      if (ncycle++ > 1000000)
         goto CLEANUP;
    
      if (tot_even)
         e = data->cycleEdgeMap[EGgetEvenPrec(curr_s,curr_t,data)];
      else
         e = data->cycleEdgeMap[EGgetOddPrec(curr_s,curr_t,data)];

      EGlistPushBack(path, e);

      e_data = (EGcycleData_t*)(e->data);

      dist = EGdijkstraCostAdd(dist, e_data->weight);

      if ( e_data->e )
         e_even = 1;
      else
         e_even = 0;

      if ( data->bdGnodeMap[e->head->id / 2]->id == curr_s )
         curr_s = data->bdGnodeMap[e->tail->id / 2]->id;
      else if ( data->bdGnodeMap[e->tail->id / 2]->id == curr_s )
         curr_s = data->bdGnodeMap[e->head->id / 2]->id;
      else if ( data->bdGnodeMap[e->head->id / 2]->id == curr_t )
         curr_t = data->bdGnodeMap[e->tail->id / 2]->id;
      else if ( data->bdGnodeMap[e->tail->id / 2]->id == curr_t )
         curr_t = data->bdGnodeMap[e->head->id / 2]->id;
      else 
      {
         fprintf(stderr, "missing case\n");
         fprintf(stderr, "curr_s = %u\n", curr_s);
         fprintf(stderr, "curr_t = %u\n", curr_t);
         fprintf(stderr, "e_head = %u\n", data->bdGnodeMap[e->head->id / 2]->id);
         fprintf(stderr, "e_tail = %u\n", data->bdGnodeMap[e->tail->id / 2]->id);
         goto CLEANUP;
      }

      if ( e_even == tot_even )
         tot_even = 1;
      else
         tot_even = 0;

   }

   if (!tot_even)
   {
      //fprintf(stderr, "need to add odd-cycle to compute actual path. IGNORING.\n");
      EGlistClear(path, nullFree);
   }
   else
   {
     EGdijkstraCost_t const evendist = EGgetEvenDistance(s,t,data);
     if ( fabs(EGdijkstraCostToLf(dist) - EGdijkstraCostToLf( evendist)) > 0.0001  )
     {
        fprintf(stderr, "computed distance  = %lf.\n", EGdijkstraCostToLf(dist));
        fprintf(stderr, "distance should be = %lf.\n", EGdijkstraCostToLf(evendist));
        fprintf(stderr, "distances do not match\n");
        rval = 1;
        goto CLEANUP;
     }
   }

   rval = 0;
   CLEANUP:

   return rval;

}

EGdkpc_t* EGdkdomToDKPC(EGdkdomino_t *dkdom,
                        EGmemPool_t *mem)
{

   EGdkpc_t *dkpc;
   unsigned int i;

   dkpc = EGmemPoolSMalloc(mem, EGdkpc_t, 1);

   dkpc->nhandles = 0;

   dkpc->dkdom = EGnewList(mem);
   for(i=0; i<EG_KPC_MAX_HANDLES; i++)
      dkpc->FH[i] = EGnewList(mem);

   EGlistPushBack(dkpc->dkdom, dkdom);

   dkpc->slack = EGdijkstraCostSub( dkdom->cut->value, EGdijkstraToCost(2.0) );

   return dkpc;

}

EGdkpc_t* EGdcutToDKPC(EGdualCut_t *dc, 
                       unsigned int orientation, 
                       EGmemPool_t *mem)
{

   EGdkdomino_t *dkdom;
   dkdom = EGdcutToDKdom( EG_KPC_MAX_HANDLES, dc, orientation, mem);
   return (EGdkdomToDKPC(dkdom, mem));

}

void EGfreeDKPC(void *v, EGmemPool_t *mem)
{

   EGdkpc_t *dkpc;
   unsigned int i;

   dkpc = (EGdkpc_t*)(v);

   for(i=0; i < EG_KPC_MAX_HANDLES; i++)
      EGfreeList(dkpc->FH[i]);

   EGlistClearMP(dkpc->dkdom, EGfreeDKdomino, mem);
   EGfreeList(dkpc->dkdom);
  
   EGmemPoolSFree(dkpc, EGdkpc_t, 1, mem);

   return;

}

EGdkdomino_t* EGdcutToDKdom(unsigned int max_k, 
                            EGdualCut_t *dc, 
                            unsigned int orientation, 
                            EGmemPool_t *mem)
{

   EGdkdomino_t *dkd=0;

   dkd = EGmemPoolSMalloc(mem, EGdkdomino_t, 1);
   dkd->k = 0;
   dkd->max_k = max_k;

   dkd->orientation=orientation;
   dkd->pairs = EGmemPoolSMalloc(mem, EGinternalPairs_t, max_k);
   dkd->handleid = EGmemPoolSMalloc(mem, unsigned int, max_k);
   dkd->cut = EGcopyDualCut(mem, dc);

   return dkd;

}

void EGfreeDKdomino(void *v, EGmemPool_t *mem)
{

   unsigned int i;
   EGdkdomino_t *dkd;

   dkd = (EGdkdomino_t*)(v);

   for(i=0; i<dkd->max_k; i++)
      EGinternalPairsClear(&(dkd->pairs[i]),mem);
   EGmemPoolSFree(dkd->pairs, EGinternalPairs_t, dkd->max_k, mem);
   EGmemPoolSFree(dkd->handleid, unsigned int, dkd->max_k, mem); 

   EGfreeDualCut(dkd->cut, mem);

   EGmemPoolSFree(dkd, EGdkdomino_t, 1, mem);
   
   return;

}

EGdijkstraCost_t EGweightDKdomino(EGdkdomino_t *dkd)
{

   unsigned int i,j;
   EGdijkstraCost_t weight;
   size_t mos[EG_MENGER_NUMOS];

   mos[EG_MENGER_ECOST]  = offsetof (EGmengerEdgeData_t, cost);
   mos[EG_MENGER_EDATA] = offsetof (EGmengerEdgeData_t, data);

   EGcomputeDualCutValue(dkd->cut);
   weight = dkd->cut->value;

   for(i=0; i < dkd->k; i++)
     for(j=0; j < dkd->pairs[i].npath; j++)
       weight = EGdijkstraCostAdd(weight, 
                EGmengerGetEdgeCost(dkd->pairs[i].stpath[j],mos));

   weight = EGdijkstraCostSub(weight, EGdijkstraToCost( (double)(dkd->k+2) ));

   return weight;

}

EGdijkstraCost_t EGcomputeDKPCslack(EGdkpc_t *dkpc)
{

   unsigned int i; 
   EGlistNode_t *it;
   size_t mos[EG_MENGER_NUMOS];

   // SLACK = 0.0

   EGdijkstraCost_t slack = EGdijkstraToCost(0.0);
   EGdijkstraCost_t weight;

   mos[EG_MENGER_ECOST]  = offsetof (EGmengerEdgeData_t, cost);
   mos[EG_MENGER_EDATA] = offsetof (EGmengerEdgeData_t, data);

   // SLACK += SUM_T w(T)

   for(it=dkpc->dkdom->begin; it; it=it->next)
   {
      weight = EGweightDKdomino( (EGdkdomino_t*)(it->this) );
      slack = EGdijkstraCostAdd(slack, weight); 
   }

   // SLACK += SUM_H x(F_H)

   for(i=0; i < dkpc->nhandles; i++)
      for(it=dkpc->FH[i]->begin; it; it=it->next)
         slack = EGdijkstraCostAdd(slack, EGmengerGetEdgeCost(it->this,mos) );

   // SLACK -= |H|
   
   slack = EGdijkstraCostSub(slack, EGdijkstraToCost( (double)(dkpc->nhandles) ) );

   return slack;

}

int EGgrowHandle(EGdkpc_t *dkpc, 
                 EGdkdomino_t *link_dkdom,
                 EGinternalPairs_t *p, 
                 EGgreedyData_t *data,
                 EGmemPool_t *mem,
                 unsigned char *const pndummy,
                 unsigned char *const pedummy)
{

   int rval, h;
   EGlist_t *extpath;
   EGlistNode_t *it;
   EGdGraphEdge_t *e;
   EGcycleData_t *edata;
   EGdkdomino_t *dkdom = 0;

   /* handle id of new handle */
   h = dkpc->nhandles;

   /* make sure we don't exceed max # of handles */
   TEST( h >= EG_KPC_MAX_HANDLES, "# of handles exceeded");

   /* extract external path */  
   extpath = p->extpath; 

   if (!extpath->size)
      goto CLEANUP;

   /* grow the dual k-domino by using the internal path */

   rval = EGgrowDKdomino(link_dkdom, p, h, mem); 
   CHECKRVAL(rval);

   /* add the handle by going through the external path */

   for(it=extpath->begin; it; it=it->next)
   {
      e = (EGdGraphEdge_t*)(it->this);
      edata = (EGcycleData_t*)(e->data);
      if (edata->e)
         EGlistPushBack(dkpc->FH[h],edata->e);
      if (edata->ddom)
      {
         dkdom = EGddominoToDKdom(edata->ddom, h, EG_KPC_MAX_HANDLES, data, mem, pndummy, pedummy);
         EGlistPushBack(dkpc->dkdom,dkdom);
      }
   }

   dkpc->nhandles += 1;
   dkpc->slack = EGdijkstraCostAdd(dkpc->slack, p->value);
   dkpc->slack = EGdijkstraCostSub(dkpc->slack, EGdijkstraToCost(2.0));

   CLEANUP:

   return 0;

}

int EGgrowDKdomino(EGdkdomino_t *dkd, 
                   EGinternalPairs_t *p,
                   int handle_num,
                   EGmemPool_t *mem)
{

   TEST( dkd->k >= dkd->max_k, "invalid dkd k size (k = %d | max = %d)", dkd->k, dkd->max_k);

   dkd->handleid[dkd->k] = handle_num;
   dkd->pairs[dkd->k].s = p->s;
   dkd->pairs[dkd->k].t = p->t;

   /* actually, the next line is misleading... 
      p->value also considers external path    */
   dkd->pairs[dkd->k].value = p->value;

   /* copy the path! */
   dkd->pairs[dkd->k].npath  = p->npath;
   dkd->pairs[dkd->k].stpath = EGmemPoolSMalloc(mem,EGdGraphEdge_t*,p->npath);
   memcpy(dkd->pairs[dkd->k].stpath,p->stpath,sizeof(EGdGraphEdge_t*)*p->npath);

   dkd->k++;

   return 0;

}

EGdkdomino_t *EGddominoToDKdom(EGddomino_t *ddom, 
                               int h, 
                               int max_nh, 
                               EGgreedyData_t *gdata,
                               EGmemPool_t *mem,
                               unsigned char *const pndummy,
                               unsigned char *const pedummy)
{
   unsigned int i;
   EGdkdomino_t *dkdom;

   dkdom = EGmemPoolSMalloc(mem, EGdkdomino_t, 1);

   dkdom->k = 1;
   dkdom->max_k = max_nh;

   dkdom->handleid = EGmemPoolSMalloc(mem,unsigned int,max_nh); 
   dkdom->handleid[0] = h;

   /* prepare internal pairs */
   dkdom->pairs = EGmemPoolSMalloc(mem,EGinternalPairs_t,max_nh); 
   dkdom->pairs[0].s      = ddom->s;
   dkdom->pairs[0].t      = ddom->t;
   dkdom->pairs[0].value  = 0;
   //dkdom->pairs[0].kdom   = dkdom;
  
   dkdom->pairs[0].npath  = ddom->npath[1];
   dkdom->pairs[0].stpath = EGmemPoolSMalloc(mem, EGdGraphEdge_t*, ddom->npath[1]);
   for(i=0; i < ddom->npath[1]; i++)
      dkdom->pairs[0].stpath[i] = ddom->path[1][i];
  
   /* dual cut */ 
   dkdom->cut = EGnewDualCut(mem, ddom->npath[0]+ddom->npath[2]);
   for(i=0; i < ddom->npath[0]; i++)
      dkdom->cut->edges[i] = ddom->path[0][i];
   for(i=0; i < ddom->npath[2]; i++)
      dkdom->cut->edges[i+ddom->npath[0]] = ddom->path[2][i];

   /* orientation */
   int rval = EGgetOrientation(gdata, 
                           dkdom->cut, 
                           dkdom->pairs[0].npath,
                           dkdom->pairs[0].stpath,
                           &(dkdom->orientation),
                           pndummy,
                           pedummy);
   if (rval)
   {
      EGfreeDKdomino(dkdom, mem);
      dkdom = 0;
   }

   return dkdom;
}

EGdualCut_t* EGextractHandleCut(EGdkpc_t *dkpc, unsigned int h)
{

   unsigned int i,j;
   EGlist_t *hcut;
   EGlistNode_t *it;
   EGdkdomino_t *dkd;
   EGinternalPairs_t *p;
   EGdualCut_t *dc;
   
   hcut = EGnewList(dkpc->dkdom->mempool); 

   for(it = dkpc->FH[h]->begin; it; it=it->next)
     EGlistPushBack(hcut, it->this);

   for(it = dkpc->dkdom->begin; it; it=it->next)
   {
      dkd = (EGdkdomino_t*)(it->this);
      for(i=0; i < dkd->k; i++)
         if ( dkd->handleid[i] == h )
         {
            p = &(dkd->pairs[i]);
            for(j=0; j<p->npath; j++)
               EGlistPushBack(hcut, p->stpath[j]); 
         }
   }

   dc = EGnewDualCut(dkpc->dkdom->mempool, hcut->size);
   for(i=0, it = hcut->begin; it; it=it->next, i++)
      dc->edges[i] = (EGdGraphEdge_t*)(it->this);
   dc->value = EGdijkstraToCost(0.0);

   EGfreeList(hcut);

   return dc;

}

void EGdisplayDualCut(FILE *fout, EGdualCut_t* dc)
{

   unsigned int i;
   size_t mos[EG_MENGER_NUMOS];

   mos[EG_MENGER_ECOST]  = offsetof (EGmengerEdgeData_t, cost);
   mos[EG_MENGER_EDATA] = offsetof (EGmengerEdgeData_t, data);

   fprintf(fout, "cut: (%d) ", dc->sz);

   for(i=0; i < dc->sz; i++)
      fprintf(fout, "(%d,%d)[%4.3lf] ", dc->edges[i]->head->id, 
              dc->edges[i]->tail->id, 
              EGdijkstraCostToLf(EGmengerGetEdgeCost(dc->edges[i],mos))); 
   fprintf(fout, "\n");

   return;

}

void EGdisplayInternalPair(FILE *fout, EGinternalPairs_t *p)
{

   unsigned int i;
   size_t mos[EG_MENGER_NUMOS];

   mos[EG_MENGER_ECOST]  = offsetof (EGmengerEdgeData_t, cost);
   mos[EG_MENGER_EDATA] = offsetof (EGmengerEdgeData_t, data);

   fprintf(fout, "(s=%d,t=%d) [sz=%d] ", p->s->id, p->t->id, p->npath);

   for(i=0; i < p->npath; i++)
      fprintf(fout, "(%d %d)[%lf] ", 
                    p->stpath[i]->tail->id, p->stpath[i]->head->id,
                    EGdijkstraCostToLf(EGmengerGetEdgeCost(p->stpath[i],mos)));

   fprintf(fout, "\n");

   return;

}

void EGdisplayDKdomino(FILE *fout, EGdkdomino_t *dkdom)
{

   unsigned int i;
   EGdijkstraCost_t weight;

   fprintf(fout, "k=%d max_k = %d [orientation=%d]\n\n", 
                 dkdom->k, dkdom->max_k, dkdom->orientation);

   fprintf(fout, "delta(T) = ");
   EGdisplayDualCut(fout, dkdom->cut);
   fprintf(fout, "\n");

   for(i=0; i<dkdom->k; i++)
   {
      fprintf(fout, "path of handle %d: ", dkdom->handleid[i]);
      EGdisplayInternalPair(fout, &(dkdom->pairs[i]));
      fprintf(fout, "\n");
   }

   weight = EGweightDKdomino(dkdom);
   fprintf(stderr, "*weight = %lf\n\n", EGdijkstraCostToLf(weight));

   return;

}


void EGdisplayDKPC(FILE *fout, EGdkpc_t *dkpc)
{

   unsigned int i,j;
   EGlistNode_t *it;
   EGdkdomino_t *dkd;
   EGdGraphEdge_t *e;

   fprintf(fout, "----------- DKPC ----------- \n\n");
   fprintf(fout, "DKPC slack    = %lf\n", EGdijkstraCostToLf(dkpc->slack));
   fprintf(fout, "nhandles = %d\n\n", dkpc->nhandles);

   for(i=0; i < dkpc->nhandles; i++)
   {
      fprintf(fout, "FH[%1u] (%3u elements): ", i, dkpc->FH[i]->size );
      for(it=dkpc->FH[i]->begin; it; it=it->next)
      {
         e = (EGdGraphEdge_t*)(it->this);
         fprintf(fout, "(%d,%d)", e->tail->id, e->head->id);
      }
      fprintf(fout, "\n");
   }

   fprintf(fout, "\n");
   fprintf(fout, "--- TEETH ---\n\n");

   for(it = dkpc->dkdom->begin; it; it=it->next)
   {
      dkd = (EGdkdomino_t *)(it->this);
      EGdisplayDKdomino(fout, dkd);
   }

   fprintf(fout, "--- HANDLES ---\n\n"); 
   
   EGdualCut_t *dc;

   for(i=0; i < dkpc->nhandles; i++)
   { 
      dc = EGextractHandleCut(dkpc, i);
      fprintf(fout, "delta(H_%d): ", i); 
      for(j=0; j < dc->sz; j++)
      {
         e = dc->edges[j];
         fprintf(fout, "(%d,%d)", e->tail->id, e->head->id);
      }
      fprintf(fout, "\n\n");
      EGfreeDualCut(dc, dkpc->dkdom->mempool);
   }
   

   return;

}

int EGprimalizeDKPC( EGdkpc_t *dkpc,
                     EGgreedyData_t *gdata,
                     EGpkpc_t *pkpc,
                     unsigned char *const pndummy,
                     unsigned char *const pedummy)
{

   int rval;
   double internal_slack,
          external_slack,
          orig_slack,
          planar_slack;
   EGdijkstraCost_t dij_external_slack;
   pkpc->randa = pkpc->randb = 0;

   rval = EGprimalizeDKPCB( dkpc,
                            gdata,
                            &pkpc->nhandles,
                            &pkpc->handle_size,
                            &pkpc->handles,
                            &pkpc->nteeth,
                            &pkpc->teeth_size,
                            &pkpc->teeth,
                            &pkpc->teeth_k,
                            &pkpc->teeth_handle,
                            &pkpc->teeth_nhalf,
                            &pkpc->teeth_halves,
                            pndummy,
                            pedummy);
   if (rval == EG_KP_NOST)
     return rval;
   CHECKRVAL(rval);

   internal_slack = (double)EGdijkstraCostToLf(dkpc->slack);
   dij_external_slack = EGcomputeDKPCslack(dkpc);
   external_slack = (double)EGdijkstraCostToLf(dij_external_slack);

   orig_slack = EGcomputePKPCslack(gdata,
                                   pkpc,
                                   gdata->norig_nodes,
                                   gdata->norig_edges,
                                   gdata->orig_edges,
                                   gdata->orig_weight);

   planar_slack = EGcomputePKPCslack(gdata,
                                     pkpc,
                                     gdata->norig_nodes,
                                     gdata->nplanar_edges,
                                     gdata->planar_edges,
                                     gdata->planar_weight);

   pkpc->slack = orig_slack;

   /* TEST FOR CONSISTENCY! */
   #if 0
   WARNING( (fabs(internal_slack - external_slack) > EG_PRIMALIZATION_ERROR) ,
         "internal slack (%lf) and external slack (%lf) are inconsistent",
          internal_slack, external_slack);
   WARNING( (fabs(internal_slack - planar_slack) > EG_PRIMALIZATION_ERROR),
         "dual slack (%lf) and primal slack (%lf) are inconsistent",
          internal_slack, planar_slack);
   #endif
   #warning A check could be added here to doubly ensure that the primal-dual conversion is correct.

   return 0;

}

int EGprimalizeDKPCB( EGdkpc_t *dkpc,
                      EGgreedyData_t *gdata,
                      int *nhandles,
                      int **handle_size,
                      int ***handles,
                      int *nteeth,
                      int **teeth_size,
                      int ***teeth,
                      int **teeth_k,
                      int ***teeth_handle,
                      int ***teeth_nhalf,
                      int ****teeth_halves,
                      unsigned char *const pndummy,
                      unsigned char *const pedummy)
{

   int rval;
   unsigned int i,j;
   EGdualCut_t *dc;
   EGlistNode_t *it;
   EGdkdomino_t *dkd;
   int kill_this_pkpc = 0;

   /* HANDLES */
   *nhandles = dkpc->nhandles;

   *handle_size = EGsMalloc(int,dkpc->nhandles);
   *handles     = EGsMalloc(int*,dkpc->nhandles);

   for(i=0; i<dkpc->nhandles; i++)
   {
      dc = EGextractHandleCut(dkpc, i);
      rval = EGgetCutNodes(gdata, 
                           dc,
                           0,
                           &((*handle_size)[i]),
                           &((*handles)[i]),
                           pndummy,
                           pedummy);
      CHECKRVAL(rval);
      EGfreeDualCut(dc, dkpc->dkdom->mempool);
      #warning Handles are not complemented before being added. 
   }

   /* TEETH */
   *nteeth = (int)(dkpc->dkdom->size); 

   *teeth_size     = EGsMalloc(int,*nteeth);
   *teeth          = EGsMalloc(int*,*nteeth);
   *teeth_k        = EGsMalloc(int,*nteeth);
   *teeth_handle   = EGsMalloc(int*,*nteeth);
   *teeth_nhalf    = EGsMalloc(int*,*nteeth);
   *teeth_halves   = EGsMalloc(int**,*nteeth);


   for(i=0, it=dkpc->dkdom->begin; it; it=it->next, i++)
   {
      dkd = (EGdkdomino_t*)(it->this);

      EGcomputeDualCutValue(dkd->cut);
     
      rval = EGgetCutNodes(gdata,
                           dkd->cut,
                           dkd->orientation,
                           &((*teeth_size)[i]),
                           &((*teeth)[i]),
                           pndummy,
                           pedummy);
      CHECKRVAL(rval);

      (*teeth_k)[i] = dkd->k;

      (*teeth_handle)[i] = EGsMalloc(int,dkd->k);
      (*teeth_nhalf)[i]  = EGsMalloc(int,dkd->k);
      (*teeth_halves)[i] = EGsMalloc(int*,dkd->k);

      for(j=0; j<dkd->k;j++)
      {
         (*teeth_handle)[i][j] = (int)dkd->handleid[j];
         rval = EGgetANodes( gdata,
                             dkd,
                             j,
                             &((*teeth_nhalf)[i][j]),
                             &((*teeth_halves)[i][j]),
                             pndummy,
                             pedummy,
                             gdata->pedges);
         if (rval == EG_KP_NOST)
            kill_this_pkpc = 1;
         else
            CHECKRVAL(rval);
      }
   }

   if (kill_this_pkpc)
      return EG_KP_NOST;

   return 0;

}

void EGfreePKPC( void *v )
{

   EGpkpc_t *pkpc = (EGpkpc_t*)(v);

   EGfreePKPCB( pkpc->nhandles,
                pkpc->handle_size,
                pkpc->handles,
                pkpc->nteeth,
                pkpc->teeth_size,
                pkpc->teeth,
                pkpc->teeth_k,
                pkpc->teeth_handle,
                pkpc->teeth_nhalf,
                pkpc->teeth_halves );

   free(pkpc);

   return;

}
                    
int EGfreePKPCB( int nhandles,
                 int *handle_size,
                 int **handles,
                 int nteeth,
                 int *teeth_size,
                 int **teeth,
                 int *teeth_k,
                 int **teeth_handle,
                 int **teeth_nhalf,
                 int ***teeth_halves )
{

   int i,j;

   for(i=0; i < nhandles; i++)
      free(handles[i]);
   
   free(handles);
   free(handle_size);

   for(i=0; i < nteeth; i++)
   {
      free(teeth[i]);    
      free(teeth_handle[i]);
      free(teeth_nhalf[i]);
      for(j=0; j < teeth_k[i]; j++)
         free(teeth_halves[i][j]);
      free(teeth_halves[i]);
   }

   free(teeth_k);
   free(teeth_size);

   free(teeth);
   free(teeth_handle);
   free(teeth_nhalf);
   free(teeth_halves);
   
   return 0;

}

int EGnewPKPCset(int nineq,
                 int **nhandles,
                 int ***handle_size,
                 int ****handles,
                 int **nteeth,
                 int ***teeth_size,
                 int ****teeth,
                 int ***teeth_k,
                 int ****teeth_handle,
                 int ****teeth_nhalf,
                 int *****teeth_halves)
{

   /* do the mallocs for nineq equations */

   *nhandles     = EGsMalloc(int,nineq);
   *handle_size  = EGsMalloc(int*,nineq);
   *handles      = EGsMalloc(int**,nineq);
   *nteeth       = EGsMalloc(int,nineq);
   *teeth_size   = EGsMalloc(int*,nineq);
   *teeth        = EGsMalloc(int**,nineq);
   *teeth_k      = EGsMalloc(int*,nineq);
   *teeth_handle = EGsMalloc(int**,nineq);
   *teeth_nhalf  = EGsMalloc(int**,nineq);
   *teeth_halves = EGsMalloc(int***,nineq);

   return 0;

}

void EGcomputeDualCutValue( EGdualCut_t *dc )
{

  unsigned int i;
  size_t mos[EG_MENGER_NUMOS];

  mos[EG_MENGER_ECOST]  = offsetof (EGmengerEdgeData_t, cost);
  mos[EG_MENGER_EDATA] = offsetof (EGmengerEdgeData_t, data);

  dc->value = EGdijkstraToCost(0.0);
  for(i=0; i < dc->sz; i++)
    dc->value = EGdijkstraCostAdd(dc->value,EGmengerGetEdgeCost(dc->edges[i],mos));

  return;

}

// dummy should be an array of size at least nnodes
double EGcomputePrimalCutValue(EGgreedyData_t*data,
                               EGpkpc_t*pkpc,
                               int set_size, 
                               int *set_nodes, 
                               int nnodes, 
                               int nedges, 
                               int *edges,
                               double *weight,
                               unsigned char *dummy)
{
   int i;
   double val = 0.0;
   memset(dummy, 0, sizeof(unsigned char)*nnodes);
   for(i=0; i < set_size; i++)
   {
      dummy[set_nodes[i]] = 1;
   }
   for(i=0; i < nedges; i++)
      if ( (dummy[edges[2*i]] + dummy[edges[2*i+1]]) == 1 )
      {
         val += weight[i];
         pkpc->randa += data->randa[i];
         pkpc->randb += data->randb[i];
      }
   return val;
}

/* dummy should be an array of size at least nnodes */

double EGcomputeInterfaceValue(int T_size, 
                               int *T_nodes, 
                               int A_size,
                               int *A_nodes,
                               int nnodes, 
                               int nedges, 
                               int *edges,
                               double *weight,
                               unsigned char *dummy)
{
   int i;
   double val = 0.0;

   memset(dummy, 0, sizeof(unsigned char)*nnodes);
  
   for(i=0; i < T_size; i++)
      dummy[T_nodes[i]] = 1;

   for(i=0; i < A_size; i++)
      if ( dummy[A_nodes[i]] )
         dummy[A_nodes[i]] = 2;
      else 
      {
         EXIT(1, "A not a subset of T");
      }

   for(i=0; i < nedges; i++)
      if ( dummy[edges[2*i]] + dummy[edges[2*i+1]] == 3 )
         val += weight[i];

   return val;
}

int EGincrementInterface( int T_size,
                          int *T,
                          int A_size,
                          int *A,
                          int nnodes,
                          int nedges,
                          int *edges,
                          unsigned char *ndummy,
                          unsigned char *emarkers )
{

   int i;
   memset(ndummy, 0, sizeof(unsigned char)*nnodes);
   for(i=0; i < T_size; i++) ndummy[ T[i] ] = 1;
   for(i=0; i < A_size; i++) ndummy[ A[i] ] = 2;
   for(i=0; i<nedges; i++)
      if ( ndummy[edges[2*i]] + ndummy[edges[2*i+1]] == 3 ) emarkers[i] += 1;
   return 0;
}

double EGcomputePKPCslack( EGgreedyData_t*data,
                           EGpkpc_t *pkpc, 
                           int nnodes,
                           int nedges,
                           int *edges,
                           double *weight )
{
  return ( EGcomputePKPCBslack( data,
                                pkpc,
                                pkpc->nhandles,
                                pkpc->handle_size,
                                pkpc->handles,
                                pkpc->nteeth,
                                pkpc->teeth_size,
                                pkpc->teeth,
                                pkpc->teeth_k,
                                pkpc->teeth_handle,
                                pkpc->teeth_nhalf,
                                pkpc->teeth_halves,
                                nnodes,
                                nedges,
                                edges,
                                weight ) );
}

double EGcomputePKPCBslack( EGgreedyData_t*data,
                            EGpkpc_t*pkpc,
                            int nhandles,
                            int *handle_size,
                            int **handles,
                            int nteeth,
                            int *teeth_size,
                            int **teeth,
                            int *teeth_k,
                            int **teeth_handle,
                            int **teeth_nhalf,
                            int ***teeth_halves,
                            int nnodes,
                            int nedges,
                            int *edges,
                            double *weight)
                           
{

   int i,j,k; 
   unsigned char *ndummy = data->pndummy, *edummy = data->pedummy;
   double slack = 0;

   /* initialize */

   /* add SUM x( delta(T) ) */
   for(i=0; i<nteeth; i++) 
      slack += EGcomputePrimalCutValue(data,
                                       pkpc,
                                       teeth_size[i],
                                       teeth[i],
                                       nnodes,
                                       nedges,
                                       edges,
                                       weight,
                                       ndummy);

   /* add SUM x(mu(H)) */
   for(i=0; i < nhandles; i++)
   {
      /* first, compute how many times each edge is in an interface */  
      memset(edummy, 0, sizeof(unsigned char)*nedges);
      for(j=0; j<nteeth;j++)
         for(k=0; k<teeth_k[j]; k++)
            if ( teeth_handle[j][k] == i )
               EGincrementInterface( teeth_size[j],
                                     teeth[j],
                                     teeth_nhalf[j][k],
                                     teeth_halves[j][k],
                                     nnodes,
                                     nedges,
                                     edges,
                                     ndummy,
                                     edummy);

      /* mark the nodes in H */
      memset(ndummy, 0, sizeof(unsigned char)*nnodes);
      for(j=0; j<handle_size[i]; j++)
         ndummy[ handles[i][j] ] = 1;

      /* add the corresponding value for each edge */
      for(j=0; j<nedges; j++)
      {
         /* if an edge is in an even number of interfaces and in delta(H) */
         if ( !(edummy[j]&1) && ( ndummy[edges[2*j]] + ndummy[edges[2*j+1]] == 1 ) )
          {
            slack += (edummy[j]+1)*weight[j];
            pkpc->randa += (edummy[j]+1)*data->randa[j];
            pkpc->randb += (edummy[j]+1)*data->randb[j];
          }

         /* if an edge is in an odd  number of interfaces and in delta(H) */
         if ( (edummy[j]&1) && ( ndummy[edges[2*j]] + ndummy[edges[2*j+1]] == 1 ) )
          {
            slack += (edummy[j])*weight[j];
            pkpc->randa += (edummy[j])*data->randa[j];
            pkpc->randb += (edummy[j])*data->randb[j];
          }

         /* if an edge 1s in an even number of interfaces and not in delta(H) */
         if ( !(edummy[j]&1) && ( ndummy[edges[2*j]] + ndummy[edges[2*j+1]] != 1 ) )
          {
            slack += (edummy[j])*weight[j];
            pkpc->randa += (edummy[j])*data->randa[j];
            pkpc->randb += (edummy[j])*data->randb[j];
          }

         /* if an edge is in an odd  number of interfaces and not in delta(H) */
         if ( (edummy[j]&1) && ( ndummy[edges[2*j]] + ndummy[edges[2*j+1]] != 1 ) )
          {
            slack += (edummy[j]+1)*weight[j];
            pkpc->randa += (edummy[j]+1)*data->randa[j];
            pkpc->randb += (edummy[j]+1)*data->randb[j];
          }
      }

   }

   /* subtract |H| */
   slack -= (double)nhandles;

   /* subtract |J| */
   for(i=0; i < nteeth; i++)
      slack -= (double)teeth_k[i];

   /* subtract 2|T| */
   slack -= (double)(2*nteeth);

   return slack;
}

int EGkpTo2pTooth( int kp_size,
                   int* kp_nodes,
                   int kp_k,
                   int* kp_handles,
                   int* kp_nhalf,
                   int** kp_halves,
                   int *tp_naset,
                   int *tp_nbset,
                   int *tp_nmset,
                   int **tp_aset,
                   int **tp_bset,
                   int **tp_mset,
                   int nnodes,
                   unsigned char *vdummy )
{

  int k, i, m[2];

  TEST( kp_k > 2 || kp_k < 1, "can only convert k-dominoes with 1 <= k <= 2");

  memset( vdummy, 0, sizeof(unsigned char)*nnodes );

  /* bit 0 == 1 iff a node is in T
     bit 1 == 1 iff a node is in A
     bit 2 == 1 iff a node is in B */

  /* mark nodes in T */
  for(i=0; i< kp_size; i++)
    vdummy[ kp_nodes[i] ] |= 1;

  /* the idea is that 
     m[k] will match handle k to handle A if kp_handle[k] = 0
     m[k] will match handle k to handle B if kp_handle[k] = 1 */

  m[0] = 2*(kp_handles[0]+1);
  m[1] = 2*(2 - kp_handles[0]);

  /* mark nodes in each half */
  for(k=0; k < kp_k; k++)
    for(i=0; i< kp_nhalf[k]; i++)
      vdummy[ kp_halves[k][i] ] |= m[k];
    
  /* count members of each set */
  *tp_naset = 0;
  *tp_nbset = 0;
  *tp_nmset = 0;
  for(i=0; i<nnodes; i++)
  {
    if (vdummy[i])
       (*tp_nmset) += 1;
    if (vdummy[i]  & 2)
       (*tp_naset) += 1;
    if (vdummy[i]  & 4)
       (*tp_nbset) += 1;
  }

  /* allocate arrays */
  if (*tp_naset)
    (*tp_aset) = EGsMalloc(int,*tp_naset);
  if (*tp_nbset)
    (*tp_bset) = EGsMalloc(int,*tp_nbset);
  if (*tp_nmset)
    (*tp_mset) = EGsMalloc(int,*tp_nmset);

  /* add members to each set */
  *tp_naset = 0;
  *tp_nbset = 0;
  *tp_nmset = 0;
  for(i=0; i<nnodes; i++)
  {
    if (vdummy[i])
       (*tp_mset)[(*tp_nmset)++] = i;
    if (vdummy[i]  & 2)
       (*tp_aset)[(*tp_naset)++] = i;
    if (vdummy[i]  & 4)
       (*tp_bset)[(*tp_nbset)++] = i;
  }

  return 0;

}

/* removes handles, very carefully, leaving the internal pairs 
   paths floating in memory - USE CAREFULLY! Not for general use */

int EGremoveHandles(EGdkpc_t *dkpc, EGmemPool_t *mem)
{

   unsigned int i;
   EGdkdomino_t *dkd;
   EGlistNode_t *it;

   /* remove the handles */
   if (!dkpc->nhandles)
      return 0;

   for(i=0; i<dkpc->nhandles; i++)
      EGlistClear(dkpc->FH[i], nullFree); 

   dkpc->nhandles = 0;

   /* redirect the paths of dkd->pairs */
   for(it=dkpc->dkdom->begin; it; it=it->next)
   {
      dkd = (EGdkdomino_t*)(it->this);
      for(i=0; i<dkd->max_k; i++)
      {
        dkd->pairs[i].npath = 0;  
        dkd->pairs[i].stpath = 0;  
      }
   }

   /* remove all the dominoes except the first one */
   while(dkpc->dkdom->size > 1)
      EGlistEraseMP(dkpc->dkdom, dkpc->dkdom->end, EGfreeDKdomino, mem);

   /* re-size the first domino */
   dkd = (EGdkdomino_t*)(dkpc->dkdom->begin->this); 
   dkd->k = 0;

   /*
   for(i=0; i<dkd->max_k; i++)
      EGinternalPairsClear(&(dkd->pairs[i]),mem);
   */

   dkpc->slack = EGcomputeDKPCslack(dkpc);

   return 0;

}

EGdualCut_t* EGcopyDualCut(EGmemPool_t *mem, EGdualCut_t *dc_in)
{

  unsigned int i;
  EGdualCut_t *dc_out;

  dc_out = EGnewDualCut(mem, dc_in->sz);

  for(i=0; i<dc_in->sz; i++)
    dc_out->edges[i] = dc_in->edges[i];

  dc_out->value = dc_in->value;

  return dc_out;

}

int KPseparateFromDualCut( EGdualCut_t *dcut,
                           unsigned int orientation,
                           EGgreedyData_t *gdata,
                           int *nadded,
                           double *best_val,
                           double sample_time )
{
 
  int rval;
  int i;
  EGmemPool_t *mem = gdata->bdG->mem;
  int add_cut_to_heap = 1;     
  double slack;

  EGdkdomino_t *zdom;
  EGdkpc_t *dkpc;
  EGpkpc_t *pkpc, *hpkpc;
  EGlist_t *pairs_list;
  EGinternalPairs_t *p, *q;
  EGlistNode_t *it, *jt;

  /* initialize */

  for( i = 0 ; i < 7 ; i++)
    if(EGdijkstraCostIsLess(dcut->value,EGdijkstraToCost(i+1.000001)))
    {
      gdata->delta_distr[i]++;
      break;
    }

  pairs_list = EGnewList(mem);

  (*nadded) = 0;
  (*best_val) = EGdijkstraToCost(1.0);
 
  /* we generate zdom only to generate the internal pairs...
     we free it immediatly after. Yes... its a bit wasteful. */
  zdom = EGdcutToDKdom( EG_KPC_MAX_HANDLES, 
                        dcut, 
                        orientation, 
                        mem);

  rval = EGgenerateInternalPairs( gdata, 
                                  zdom, 
                                  zdom->orientation, 
                                  pairs_list, 
                                  (unsigned int)gdata->max_handles*15, 
                                  gdata->percentage,
                                  gdata->pndummy,
                                  gdata->pedummy,
                                  gdata->bdnodes);
  CHECKRVAL(rval);

  EGfreeDKdomino(zdom, mem);
  zdom = 0;

  if (sample_time > 0.1)
  {
    //fprintf(stdout, "before: %u  ", pairs_list->size);
    EGlist_t *new_pairs;
    // generate the new list
    new_pairs = EGnewList(gdata->bdG->mem);
    rval = EGgenerateRandomInternalPairs( gdata,
                                          new_pairs,
                                          pairs_list,
                                          (unsigned int)gdata->max_handles*35, 
                                          sample_time );
    CHECKRVAL(rval);
    // merge the lists
    pairs_list = EGmergeIPlists(pairs_list, new_pairs);
    //fprintf(stdout, "after: %u      ", pairs_list->size);
  }

  /* now we loop through all possible combinations of internal
     paths so as to generate as many constraints as possible   */
  for(it=pairs_list->begin; it; it=it->next)
  for(jt=it->next; jt; jt=jt->next)
  {
        
    p = (EGinternalPairs_t*)(it->this);
    q = (EGinternalPairs_t*)(jt->this);

    if (p->s == q->s && p->t == q->t)
      continue;

    zdom = EGdcutToDKdom( EG_KPC_MAX_HANDLES, 
                          dcut, 
                          orientation, 
                          mem);

    dkpc = EGdkdomToDKPC(zdom, mem);

    slack = EGdijkstraCostToLf(dkpc->slack);


    rval = EGgrowHandle(dkpc, 
                        zdom,
                        p, 
                        gdata, 
                        mem,
                        gdata->pndummy,
                        gdata->pedummy);
    CHECKRVAL(rval); 

    rval = EGgrowHandle(dkpc, 
                        zdom,
                        q, 
                        gdata, 
                        mem,
                        gdata->pndummy,
                        gdata->pedummy);
    CHECKRVAL(rval);

    /* if the current constraint is no good, then, because
    the pairs are sorted, the next ones might not be as well */
    if (EGdijkstraCostToLf(dkpc->slack) > -0.01)
    {
      EGfreeDKPC(dkpc, mem);
      dkpc = 0;
      if (jt == pairs_list->begin->next)
        it = pairs_list->end;
      break;
    }

    /* primalize the constraint */
    pkpc = EGsMalloc(EGpkpc_t,1);

    add_cut_to_heap = 1;

    rval = EGprimalizeDKPC( dkpc,
                            gdata,
                            pkpc,
                            gdata->pndummy,
                            gdata->pedummy);
    if (rval == EG_KP_NOST)
    {
       rval = 0;
       add_cut_to_heap = 0;
    }
    else if (rval)
    {
      fprintf(stdout, "WARNING! graphs are inconsistent. Saving to 'inconsistent_graph.b.x'\n");
      saveBgraph ("inconsistent_graph.b.x", gdata->norig_nodes, 
                                            gdata->norig_edges, 
                                            gdata->orig_edges, 
                                            gdata->orig_weight);
    }
    CHECKRVAL(rval);

    /* check if the cut is violated */
   if (pkpc->slack > 0.0001)
     add_cut_to_heap = 0;

    /* check if handle is empty */
    for(i=0; i<pkpc->nhandles; i++)
      if ( !pkpc->handle_size[i] || (pkpc->handle_size[i] == gdata->norig_nodes) )
      {
        #warning Erasing constraints with empty handles.
        //MESSAGE(0, "Warning! Removing a constraint with empty handle!");
        add_cut_to_heap = 0;
        break;
      }
    
    /* if the heap is full, check if the cut has fabs(slack) 
    lower than all the elements of the heap               */
    if (gdata->cut_heap->size == gdata->cut_heap->max_size)
    {
      const EGheapCost_t mv = EGheapGetMinVal(gdata->cut_heap);
      const double mv_lf = (double)EGheapCostToLf(mv);
      if ( !(fabs(pkpc->slack) > mv_lf) )
        add_cut_to_heap = 0;
    }
 
    /* check if the cut is already in the constraint heap */
    if (add_cut_to_heap)
      for(i=0; i < (int)gdata->cut_heap->size; i++)
      {
        hpkpc = (EGpkpc_t*)(gdata->cut_heap->heap[i]->this);
        if (hpkpc->randa == pkpc->randa)
          if (hpkpc->randb == pkpc->randb) 
            if (hpkpc->nhandles == pkpc->nhandles)
              if (hpkpc->nteeth == pkpc->nteeth)
                add_cut_to_heap = 0;   
      } 

    if (add_cut_to_heap)
    {
      /* if the cut heap is full, remove constraint with highest slack */
      if (gdata->cut_heap->size == gdata->cut_heap->max_size)
      {
        hpkpc = (EGpkpc_t*)(EGheapGetMin(gdata->cut_heap)->this);
        EGfreePKPC(hpkpc);
        EGheapDeleteMin(gdata->cut_heap, mem);
      }

      /* add the new constraint to the heap */
      EGheapInsert( gdata->cut_heap,
                    pkpc,
                    EGheapToCost( fabs(pkpc->slack) ),
                    mem );
      if ( !(*nadded) )
         (*best_val) = pkpc->slack;
      (*nadded) += 1;
    } 


    /* update global statistics */
    if (pkpc->slack < -0.01)
      gdata->nviolated += 1;
    if (pkpc->slack < gdata->min_violated)
        gdata->min_violated = pkpc->slack;
    if (pkpc->slack > gdata->max_violated)
        gdata->max_violated = pkpc->slack;

    #if DISPLAY_MODE
    if (add_cut_to_heap)
      EGdisplayStatus(stdout, gdata);
    #endif

    /* free dual form */
    EGfreeDKPC(dkpc, mem);
    dkpc = 0;

    /* if we do not keep the primal cut, we free it */
    if (!add_cut_to_heap && pkpc)
      EGfreePKPC(pkpc);
    pkpc = 0;

  }

  /* free objects used to derive this round of cuts */
  EGlistClearMP(pairs_list,EGfreeInternalPairs,mem);
  EGfreeList(pairs_list);

  return 0;

}
                           
void EGdisplayStatus(FILE *fout, EGgreedyData_t *gdata)
{

  double hmax;
  EGheapCost_t hhmax;

  if (gdata->cut_heap->size)
  {
    hhmax = EGheapGetMinVal(gdata->cut_heap);
    hmax = EGheapCostToLf(hhmax);
  }
  else
    hmax = 0.0;

  //if(!(gdata->nviolated&0x1f))
  if(!(gdata->nviolated&0xf))
  fprintf(fout, "\rnviolated = %4u   min = %10lf   max = %10lf   hsz = %3u   hmax = %10lf", 
          gdata->nviolated, gdata->min_violated, gdata->max_violated, 
          gdata->cut_heap->size, -1*hmax);

  return;

}


int KPprocess_cut(int*cutset,
                  int cutsize,
                  double cutweight,
                  void*process_indo,
                  setlist**sets)
{
  int rval = 0, nadded = 0;
  EGgreedyData_t*const data = (EGgreedyData_t*) process_indo;
  EGdualCut_t*dualcut = 0;
  double minslack;
  EGmemPool_t *mem = data->bdG->mem;
  
  /* to eliminate warnings from the compiler */
  sets = 0;
  cutweight = 0;
  
  /* get the dual cut */
  rval = EGgetDualCut(data, &dualcut, cutsize, cutset, data->pndummy,
                      data->pedummy, data->pedges);
  CHECKRVALG(rval,CLEANUP);

  /* separate: orientation = 0 */
  rval = KPseparateFromDualCut(dualcut, 0, data, &nadded, &minslack, 0.0);
  CHECKRVALG(rval,CLEANUP);

  if (nadded)
  {
    EGdualCut_t *dc_copy = EGcopyDualCut(mem, dualcut);
    EGaddDualCutToHeap(data->dc_heap,
                       dc_copy,
                       0, 
                       minslack,
                       mem);
  }

  /* separate: orientation = 1 */
  rval = KPseparateFromDualCut(dualcut, 1, data, &nadded, &minslack, 0.0);
  CHECKRVALG(rval,CLEANUP);

  if (nadded)
  {
    EGdualCut_t *dc_copy = EGcopyDualCut(mem, dualcut);
    EGaddDualCutToHeap(data->dc_heap,
                       dc_copy,
                       1, 
                       minslack,
                       mem);
  }

  EGfreeDualCut(dualcut,data->bdG->mem);

  /* ending */
  CLEANUP:
  return rval ? PROCESS_ERROR:PROCESS_INCOMPLETE;
}

int EGaddDualCutToHeap(EGheap_t *h, 
                       EGdualCut_t *dc, 
                       int orientation,
                       double lfval, 
                       EGmemPool_t *mem)
{

  unsigned int i;
  EGheapCost_t hval = EGheapToCost(fabs(lfval));
  EGdualCut_t *hdc;
  EGcutSeed_t *hcs;
  int add_dc_to_heap = 0;

  /* sort dual cut data */
  qsort (dc->edges, dc->sz, sizeof (EGdGraphEdge_t *), EGedgeCompare);

  /* check first if the dual cut is already in the heap */
  for(i=0; i<h->size; i++)
  {
    hcs = (EGcutSeed_t*) h->heap[i]->this;
    hdc = hcs->dc;
    if ( hcs->orientation == orientation && EGareDualCutsEqual(hdc,dc) )
      goto DONE; 
  }

  /* if the heap is not full, in goes the cut */
  if ( h->size != h->max_size )
    add_dc_to_heap = 1;
  /* else, we compare it to the worst value in the heap */
  else if (EGheapCostIsLess(EGheapGetMinVal(h),hval))
  {
    /* if its better, we remove this 'bad' heap entry */
    hcs = (EGcutSeed_t*)EGheapGetMinThis(h);
    EGfreeDualCut(hcs->dc,mem);
    EGmemPoolSFree(hcs, EGcutSeed_t, 1, mem);
    EGheapDeleteMin(h,mem);
    add_dc_to_heap = 1;
  }

  DONE:

  /* add to the heap if we have decided to do so */
  if (add_dc_to_heap)
  {
    EGcutSeed_t *cs;
    cs = EGmemPoolSMalloc(mem, EGcutSeed_t, 1);
    cs->dc = dc;
    cs->orientation = orientation;
    EGheapInsert(h,cs,hval,mem);
  }
  else
    EGfreeDualCut(dc,mem);
  
  return 0; 

}

int EGareDualCutsEqual(EGdualCut_t *a, EGdualCut_t *b)
{
  unsigned int i, are_equal = 0;

  if (a->sz != b->sz)
    goto DONE;

  for (i = 0; i < a->sz; i++)
    if (EGboyerEdgeNumber (a->edges[i]) != EGboyerEdgeNumber (b->edges[i]))
      goto DONE;

  are_equal = 1;

  DONE:

  return are_equal;
}

---