eg_kppairs.c

Go to the documentation of this file.

#include "eg_kppairs.h"

#define GETOTHEREND(_cedge) (G->G[gp_GetTwinArc(G,_cedge - G->G)].v)
#define GETROOTEDGE(curroot) (G->G + G->G[curroot].link[1])
#define GETNEXTEDGE(curedge,curroot) ((curedge->link[1] < G->N) ? GETROOTEDGE(curroot) : (G->G + curedge->link[1]))

/* ========================================================================= */
#define EG_KPPAIRS_VERBOSE 100
#define EG_KPPAIRS_DEBUG 0

static size_t os[EG_MENGER_NUMOS] = {
    [EG_MENGER_PI] = offsetof (EGmengerNodeData_t, pi),
    [EG_MENGER_DIST] = offsetof (EGmengerNodeData_t, dist),
    [EG_MENGER_ORIG_DIST] = offsetof (EGmengerNodeData_t, orig_dist),
    [EG_MENGER_NDIST] = offsetof (EGmengerNodeData_t, ndist),
    [EG_MENGER_ORIG_NDIST] = offsetof (EGmengerNodeData_t, orig_ndist),
    [EG_MENGER_MARK] = offsetof (EGmengerNodeData_t, marker),
    [EG_MENGER_ORIG_MARK] = offsetof (EGmengerNodeData_t, orig_marker),
    [EG_MENGER_FATHER] = offsetof (EGmengerNodeData_t, father),
    [EG_MENGER_ORIG_FATHER] = offsetof (EGmengerNodeData_t, orig_father),
    [EG_MENGER_HEAP_CONNECTOR] = offsetof (EGmengerNodeData_t, hc),
    [EG_MENGER_ECOST] = offsetof (EGmengerEdgeData_t, cost),
    [EG_MENGER_REDUCED_ECOST] = offsetof (EGmengerEdgeData_t, reduced_cost),
    [EG_MENGER_IS_IN_SOL] = offsetof (EGmengerEdgeData_t, is_in_solution),
    [EG_MENGER_EDATA] = offsetof (EGmengerEdgeData_t, data)
  };
static size_t  dij_os[6] = { 
    [EG_DIJ_DIST] = offsetof (EGmengerNodeData_t, orig_dist),
    [EG_DIJ_NDIST] = offsetof (EGmengerNodeData_t, orig_ndist),
    [EG_DIJ_FATHER] = offsetof (EGmengerNodeData_t, orig_father),
    [EG_DIJ_MARKER] = offsetof (EGmengerNodeData_t, orig_marker),
    [EG_DIJ_HCONNECTOR] = offsetof (EGmengerNodeData_t, hc),
    [EG_DIJ_ELENGTH] = offsetof (EGmengerEdgeData_t, cost)
  };

#define KP_CHECK_SUM_IS_LESS(a,b,c) \
  WARNINGL(EG_KPPAIRS_DEBUG, EGdijkstraCostIsLess( EGdijkstraCostAdd((a),(b)), EGdijkstraToCost(c)), "Inequality not satisfied %lf < %lf", EGdijkstraCostToLf(EGdijkstraCostAdd((a), (b))), c)
/* ========================================================================= */
int static EGdoPrimalDfs( const unsigned n_marks,
                          size_t*const set_sz,
                          EGlist_t** CUT, 
                          unsigned char*const node_mark, 
                          unsigned char*const edge_mark, 
                          EGgreedyData_t*const data)
{
  int rval = 0;
  const size_t nedges = data->bdG->edge_id_max/2;
  const size_t nnodes = nedges + 2 - data->bdG->node_id_max;
  size_t node_id = 0, other_id;
  const graphP G = data->G;
  register unsigned i;
  unsigned char cur_mark = 0, oth_mark;
  graphNodeP pedge;
  for( i = n_marks ; i-- ; )
  {
    set_sz[i] = 0;
    EGlistClear(CUT[i],0);
  }
  memset(node_mark,0,sizeof(unsigned char)*nnodes);
  node_mark[node_id] = cur_mark+1;
  set_sz[cur_mark] += 1;
  EGlistPushHead(CUT[cur_mark],(void*)node_id);
  do
  {
    node_id = (size_t)(CUT[cur_mark]->begin->this);
    EGlistErase(CUT[cur_mark],CUT[cur_mark]->begin,0);
    pedge = GETROOTEDGE(node_id);
    do
    {
      if(edge_mark[pedge->id]&128U) goto NEXT_EDGE;
      edge_mark[pedge->id] |= 128U;
      other_id = pedge->v;
      oth_mark = cur_mark ^ ((edge_mark[pedge->id])&127U);
      if(node_mark[other_id])
      {
        TESTG((rval=(oth_mark+1!=node_mark[other_id])), CLEANUP, "Imposible");
        goto NEXT_EDGE;
      }
      node_mark[other_id] = oth_mark + 1;
      set_sz[oth_mark] += 1;
      EGlistPushHead(CUT[oth_mark],(void*)other_id);
      NEXT_EDGE:
      pedge = GETNEXTEDGE(pedge, node_id);
    } while(pedge != GETROOTEDGE(node_id));
    cur_mark = 0;
    while( cur_mark < n_marks && !(CUT[cur_mark]->size)) cur_mark++;
  } while( cur_mark < n_marks);
  node_id = 0;
  #if KPPAIRS_VERBOSE <= DEBUG
  fprintf(stderr,"Number of Sets found %u\n",n_marks);
  #endif
  for( i = n_marks ; i-- ; )
  {
    node_id += set_sz[i];
    #if KPPAIRS_VERBOSE <= DEBUG
    fprintf(stderr,"\tset[%u]=%zd\n",i,set_sz[i]);
    #endif
  }
  TESTG((rval=(node_id != nnodes)), CLEANUP, "Imposible, nnodes %zd, marked nodes %zd",nnodes,node_id);
  CLEANUP:
  return rval;
}


/* ========================================================================= */
int EGgenerateInternalPairs (EGgreedyData_t*const data,
                             EGdkdomino_t*const zerodom,
                             const unsigned int orientation,
                             EGlist_t*const pairs,
                             const unsigned max_pairs,
                             double const percentage,
                             unsigned char*const node_mark,
                             unsigned char*const edge_mark,
                             EGdGraphNode_t**const dc_nodes)
{
  unsigned register int i = 0, j = 0, k = 0;
  int node_id, other_id;
  //EGlist_t**dembed = data->dembed;
  EGdGraph_t*const bdG = data->bdG;
  EGmemPool_t* mem = bdG->mem;
  EGlist_t*CUT[2]={EGnewList(mem),EGnewList(mem)};
  size_t set_sz[2] = {0,0};
  EGinternalPairs_t *worst_pair = 0, *cur_pair = 0;
  unsigned char const cur_mark = orientation ? 2U : 1U;
  unsigned int tot_pairs = 0;
  EGlist_t* unattach_edges = EGnewList(mem);
  EGdGraphNode_t *cur_node;
  EGdGraphEdge_t *tmp_edge = 0;
  EGheap_t* dij_heap = EGnewHeap(mem, 2, bdG->nodes->size);
  EGlistNode_t* It1;
  const graphP G = data->G;
  graphNodeP pedge;
  unsigned int cycle_sz = zerodom->cut->sz;
  int rval = 0;
  EGlist_t *extpath;
  MESSAGE(KPPAIRS_VERBOSE,"Entering for cut %p orientation %u", 
          (void*)zerodom->cut, orientation);
  
  EGdijkstraCost_t dist, path_dist;
  EGdijkstraCost_t const slack = EGdijkstraCostMul(EGdijkstraToCost(2.0),
                                                EGdijkstraToCost(percentage));
  /* test that the graph max id is equal to the number of nodes in it */
  TESTG((rval=(bdG->node_id_max!=bdG->nodes->size)),CLEANUP,"Imposible!");

  /* now we define the cut as the edges with mark 1. and identify all nodes 
   * in the (dual) cycle */
  memset(dc_nodes,0,sizeof(EGdGraphNode_t*)*bdG->nodes->size);
  memset(edge_mark,0,sizeof(unsigned char)*bdG->edge_id_max/2);
  for( i = cycle_sz ; i-- ; )
  {
    WARNINGL(KPPAIRS_VERBOSE, edge_mark[ EGmengerGetEdgeData( 
             zerodom->cut->edges[i],os)->id], "repeated edge %p in cut", 
             (void*)zerodom->cut->edges[i]);
    edge_mark[EGmengerGetEdgeData(zerodom->cut->edges[i],os)->id] ^= 1U;
    dc_nodes[zerodom->cut->edges[i]->head->id] = zerodom->cut->edges[i]->head;
    dc_nodes[zerodom->cut->edges[i]->tail->id] = zerodom->cut->edges[i]->tail;
  }

  /* identify the primal partition */
  rval = EGdoPrimalDfs(2,set_sz,CUT,node_mark,edge_mark,data);
  CHECKRVALG(rval,CLEANUP);
  TESTG((rval=!set_sz[cur_mark-1]),CLEANUP,"Empty set! imposible!");

  /* put all nodes in the delta of the cycle at the beginning of the array */
  for( i = 0, tot_pairs = 0 ; i < bdG->nodes->size ; i++ )
    if(dc_nodes[i]) dc_nodes[tot_pairs++] = dc_nodes[i];
  TESTG((rval=(cycle_sz < tot_pairs)),CLEANUP,"Imposible!, cycle_sz %u "
        "cycle_nodes %u",cycle_sz, tot_pairs);
  cycle_sz = tot_pairs;

  /* depending in the orientation, we unnatach the edges with at least one
   * primal end in the wrong side for all nodes in the boundary of the cycle. 
   * */
  for( i = cycle_sz ; i-- ; )
  {
    cur_node = dc_nodes[i];
    MESSAGE(KPPAIRS_VERBOSE,"Checking node %u", cur_node->id);
    for( It1 = cur_node->out_edges->begin ; It1 ; It1 = It1->next)
    {
      pedge = EGmengerGetEdgeData(It1->this,os);
      node_id = pedge->v;
      other_id = GETOTHEREND(pedge);
      if(node_mark[node_id] != cur_mark || node_mark[other_id] != cur_mark)
      {
        MESSAGE(KPPAIRS_VERBOSE,"Unnattaching edge %p: %u (%u,%u)", 
                It1->this, ((EGdGraphEdge_t*)It1->this)->id, 
                ((EGdGraphEdge_t*)It1->this)->tail->id, 
                ((EGdGraphEdge_t*)It1->this)->head->id);
        EGlistPushBack(unattach_edges,It1->this);
        //EGdGraphUnattachEdge(bdG,((EGdGraphEdge_t*)It1->this));
      }
    }
  }
  for(It1 = unattach_edges->begin ; It1; It1 = It1->next)
    EGdGraphUnattachEdge(bdG,((EGdGraphEdge_t*)It1->this));
  /* display the reduced graph */
  #if KPPAIRS_VERBOSE <= DEBUG
  MESSAGE(KPPAIRS_VERBOSE,"Unattach %u edges from G", unattach_edges->size);
  EGdGraphDisplay(bdG,0,0,0,stderr);
  #endif

  /* now we are ready to call dijkstra in the reduced graph */
  tot_pairs = 0;
  for( i = cycle_sz ; i-- ; )
  {
    rval = EGpartialDijkstra( dc_nodes[i], 0, slack , dij_os, dij_heap, bdG);
    CHECKRVALG(rval,CLEANUP);
    for( j = i ; j-- ; )
    {
      /* discard pairs that are unreachable */
      if(EGmengerGetOrigMark(dc_nodes[j],os) > UINT_MAX - 2) continue;
      /* discard pairs with cost over the bound */
      path_dist = EGdijkstraGetDist(dc_nodes[j],dij_os);
      dist = EGdijkstraCostAdd( path_dist, EGgetEvenDistance(dc_nodes[i]->id, 
                                dc_nodes[j]->id, data));
      /* note that cut->val + dist < 4 for the inequality to have a chance to
       * be violated, se we ask that dist < 4 - cut->val <==> dist < slack,
       * thus we discard all pairs such that slack < dist */
      if(EGdijkstraCostIsLess(slack,dist)) continue;
      /* discard pairs that are already in the kdomino */
      for(k = zerodom->k ; k-- ; )
      {
        if(dc_nodes[i] == zerodom->pairs[k].s && 
           dc_nodes[j] == zerodom->pairs[k].t) 
          continue;
        if(dc_nodes[i] == zerodom->pairs[k].t && 
           dc_nodes[j] == zerodom->pairs[k].s) 
          continue;
      }
      MESSAGE(KPPAIRS_VERBOSE,"Path (%u,%u) with value %lf", dc_nodes[i]->id,
              dc_nodes[j]->id, EGdijkstraCostToLf(dist));
      /* compare against worst pair */
      if(pairs->size >= max_pairs)
      {
        worst_pair = (EGinternalPairs_t*)pairs->end->this;
        if(EGdijkstraCostIsLess(worst_pair->value,dist)) continue;
      }

      /* ensure we can extract the even path */
      extpath = EGnewList(mem);
      EGextractEEpath(dc_nodes[i]->id, dc_nodes[j]->id, data, extpath);
      if (!extpath->size)
      {
        #warning WARNING: This happens a lot!! Look.
        EGlistClear(extpath, nullFree);
        EGfreeList(extpath);
        continue;
      }

      /* define the current pair */
      tot_pairs++;
      cur_pair = EGmemPoolSMalloc(mem,EGinternalPairs_t,1);
      cur_pair->s = dc_nodes[i];
      cur_pair->t = dc_nodes[j];
      cur_pair->ivalue = path_dist;
      cur_pair->value = dist;
      cur_pair->stpath = EGmemPoolSMalloc(mem, EGdGraphEdge_t*, EGdijkstraGetNdist(dc_nodes[j], dij_os));
      EGdijkstraExtractSolution( cur_pair->stpath, &cur_pair->npath, dc_nodes[i], dc_nodes[j], dij_os );
      cur_pair->extpath = extpath;

      /* otherwise, we add the pair to the external heap of pairs */
      for(It1 = pairs->end ; It1 ; It1 = It1->prev)
      {
        worst_pair = (EGinternalPairs_t*)It1->this;
        if(EGdijkstraCostIsLess(worst_pair->value,dist)) 
          break;
      }
      if(It1) EGlistInsertAfter(pairs,It1,cur_pair);
      else EGlistPushHead(pairs,cur_pair);
      while(pairs->size > max_pairs)
      {
        EGlistEraseMP(pairs,pairs->end,EGfreeInternalPairs,mem);
      }
    }
  }
  
  /* ending */
  CLEANUP:
  if(tot_pairs)
  MESSAGE(KPPAIRS_VERBOSE, "for kdom %p generate %u pairs,  value range"
          " [%lf,%lf], eliminated edges %u",
          (void*)zerodom, tot_pairs, pairs->size ? 
          EGdijkstraCostToLf(((EGinternalPairs_t*)pairs->begin->this)->value):
          2.0, pairs->size ? 
          EGdijkstraCostToLf(((EGinternalPairs_t*)pairs->end->this)->value):
          2.0, unattach_edges->size);
  for(It1 = unattach_edges->begin ; It1; It1 = It1->next)
  {
    tmp_edge = (EGdGraphEdge_t*)It1->this;
    EGdGraphAttachEdge(bdG,tmp_edge);
  }
  EGfreeList(unattach_edges);
  EGfreeHeap(dij_heap,mem);
  EGfreeList(CUT[0]);
  EGfreeList(CUT[1]);
  return rval;
}

/* ========================================================================= */
void EGinternalPairsClear(void*pair,EGmemPool_t*mem)
{
 
  EGinternalPairs_t*const p = (EGinternalPairs_t*)pair;
  if(p->stpath) EGmemPoolSFree(p->stpath, EGdGraphEdge_t*, p->npath, mem); 
  p->stpath = 0;
  p->npath = 0;
  if(p->extpath)
  {
    EGlistClear(p->extpath, nullFree);
    EGfreeList(p->extpath);
    p->extpath = 0;
  }
  return;
}

/* ========================================================================= */
void EGfreeInternalPairs(void*pair,EGmemPool_t*mem)
{
  EGinternalPairsClear(pair,mem);
  EGmemPoolSFree(pair,EGinternalPairs_t,1,mem);
  return;
}

/* ========================================================================= */
int EGgetCutNodes(EGgreedyData_t*const data,
                  EGdualCut_t*const dualcut,
                  const unsigned int orientation,
                  int*const cutsz,
                  int**const cutnodes,
                  unsigned char*const node_mark,
                  unsigned char*const edge_mark)
{
  unsigned register int i = 0;
  int j = 0;
  EGmemPool_t* mem = data->bdG->mem;
  EGlist_t*CUT[2]={EGnewList(mem),EGnewList(mem)};
  size_t set_sz[2] = {0,0};
  unsigned char cur_mark = orientation ? 2U : 1U;
  unsigned const int cycle_sz = dualcut->sz;
  const size_t nedges = data->bdG->edge_id_max/2;
  const size_t nnodes = nedges + 2 - data->bdG->node_id_max;
  int rval = 0;
  
  /* test that the graph max id is equal to the number of nodes in it */
  TESTG((rval=(data->bdG->node_id_max!=data->bdG->nodes->size)), 
        CLEANUP,"Imposible!");

  /* now we define the cut as the edges with mark 1. and identify all nodes 
   * in the (dual) cycle */
  memset(edge_mark,0,sizeof(unsigned char)*nedges);
  for( i = cycle_sz ; i-- ; )
    edge_mark[EGmengerGetEdgeData(dualcut->edges[i],os)->id] ^= 1U;

  /* identify the partition */
  rval = EGdoPrimalDfs(2,set_sz,CUT,node_mark,edge_mark,data);
  CHECKRVALG(rval,CLEANUP);
  TESTG((rval=!set_sz[cur_mark-1]),CLEANUP,"Empty set! imposible!");

  /* depending on the orientation, we select the size of the set */
  *cutsz = set_sz[cur_mark-1];
  if(*cutsz) *cutnodes = EGsMalloc(int,*cutsz);
  j = 0;
  for( i = nnodes ; i-- ; )
    if(node_mark[i] == cur_mark)
      (*cutnodes)[j++] = i;
  TESTG((rval=(j!=(*cutsz))), CLEANUP, "Imposible!");

  /* ending */
  CLEANUP:
  EGfreeList(CUT[0]);
  EGfreeList(CUT[1]);
  return rval;
}

/* ========================================================================= */
int EGgetANodes(EGgreedyData_t*const data,
                EGdkdomino_t*const kdom,
                int unsigned const k,
                int*const cutsz,
                int**const cutnodes,
                unsigned char*const node_mark,
                unsigned char*const edge_mark,
                EGdGraphEdge_t**const dc_nodes)
{
  unsigned register int i = 0;
  int j = 0;
  EGmemPool_t* mem = data->bdG->mem;
  EGlist_t*CUT[4]={EGnewList(mem),EGnewList(mem),EGnewList(mem),EGnewList(mem)};
  size_t set_sz[4] = {0,0,0,0};
  unsigned char cur_mark = kdom->orientation ? 4U : 2U;
  unsigned const int cycle_sz = kdom->cut->sz;
  const size_t nedges = data->bdG->edge_id_max/2;
  const size_t nnodes = nedges + 2 - data->bdG->node_id_max;
  EGdGraphNode_t*s = kdom->pairs[k].s;
  EGdGraphNode_t*const t = kdom->pairs[k].t;
  EGdGraphEdge_t*e=0;
  int rval = 0;
  EGlistNode_t* It = 0;
  *cutsz = 0;
  *cutnodes = 0;
  EGlist_t* DFS = EGnewList(mem);
  memset(dc_nodes,0,sizeof(EGdGraphEdge_t*)*data->bdG->nodes->size);
  
  /* test that the graph max id is equal to the number of nodes in it */
  TESTG((rval=(data->bdG->node_id_max!=data->bdG->nodes->size)), 
        CLEANUP,"Imposible!");

  /* now we define the cut as the edges with mark 2. and identify all nodes 
   * in the (dual) cycle */
  memset(edge_mark,0,sizeof(unsigned char)*nedges);
  for( i = cycle_sz ; i-- ; )
    edge_mark[EGmengerGetEdgeData(kdom->cut->edges[i],os)->id] ^= 2U;

  /* now we must find another s-t path with mark 2 and set the mark ^= 1U */
  EGlistPushBack(DFS,s);
  while(DFS->size && s!= t)
  {
    s = (EGdGraphNode_t*)DFS->begin->this;
    EGlistErase(DFS,DFS->begin,0);
    for( It = s->out_edges->begin ; It ; It = It->next)
    {
      e = (EGdGraphEdge_t*)It->this;
      if(edge_mark[EGmengerGetEdgeData(e,os)->id]== 2 && !dc_nodes[e->head->id])
      {
        dc_nodes[e->head->id] = e;
        EGlistPushHead(DFS,e->head);
        edge_mark[EGmengerGetEdgeData(e,os)->id] = 3;
      }
    }
  }
  if(s!=t)
  {
    rval = EG_KP_NOST;
    //MESSAGE(0, "Imposible, can't reach T from S");
    goto CLEANUP;
  }

  /* now we define the cut as the edges with mark 2. and identify all nodes 
   * in the (dual) cycle */
  memset(edge_mark,0,sizeof(unsigned char)*nedges);
  for( i = cycle_sz ; i-- ; )
    edge_mark[EGmengerGetEdgeData(kdom->cut->edges[i],os)->id] ^= 2U;

  /* now we mark all edges in the path with mark 1 */
  for( i = kdom->pairs[k].npath ; i-- ; )
    edge_mark[EGmengerGetEdgeData(kdom->pairs[k].stpath[i],os)->id] ^= 1U;
  
  /* now we mark the other path from s to t */
  for( e = dc_nodes[t->id] ; e && e->head != kdom->pairs[k].s ; e = dc_nodes[e->tail->id])
    edge_mark[EGmengerGetEdgeData(e,os)->id] ^= 1U;

  /* identify the partition */
  rval = EGdoPrimalDfs(4,set_sz,CUT,node_mark,edge_mark,data);
  CHECKRVALG(rval,CLEANUP);
  TESTG((rval=!set_sz[cur_mark-1]),CLEANUP,"Empty set! imposible!");

  /* depending on the orientation, we select the size of the set */
  *cutsz = set_sz[cur_mark-1];
  if(*cutsz) *cutnodes = EGsMalloc(int,*cutsz);
  j = 0;
  for( i = nnodes ; i-- ; )
    if(node_mark[i] == cur_mark)
      (*cutnodes)[j++] = i;
  TESTG((rval=(j!=(*cutsz))), CLEANUP, "Imposible!");

  /* ending */
  CLEANUP:
  EGfreeList(DFS);
  EGfreeList(CUT[0]);
  EGfreeList(CUT[1]);
  EGfreeList(CUT[2]);
  EGfreeList(CUT[3]);
  return rval;
}

/* ========================================================================= */
int EGgetOrientation( EGgreedyData_t*const data,
                      EGdualCut_t*const dualcut,
                      int unsigned const pathsz,
                      EGdGraphEdge_t**const path,
                      unsigned int*const orientation,
                      unsigned char*const node_mark,
                      unsigned char*const edge_mark)
{
  unsigned register int i = 0;
  EGmemPool_t* mem = data->bdG->mem;
  EGlist_t*CUT[2]={EGnewList(mem),EGnewList(mem)};
  size_t set_sz[2] = {0,0};
  unsigned char cur_mark = 0;
  unsigned const int cycle_sz = dualcut->sz;
  const size_t nedges = data->bdG->edge_id_max/2;
  const graphP G = data->G;
  graphNodeP pedge;
  int rval = 0;
  
  /* test that the graph max id is equal to the number of nodes in it */
  TESTG((rval=(data->bdG->node_id_max!=data->bdG->nodes->size)), 
        CLEANUP,"Imposible!");

  /* now we define the cut as the edges with mark 1. and identify all nodes 
   * in the (dual) cycle */
  memset(edge_mark,0,sizeof(unsigned char)*nedges);
  for( i = cycle_sz ; i-- ; )
    edge_mark[EGmengerGetEdgeData(dualcut->edges[i],os)->id] ^= 1U;

  /* identify the partition */
  rval = EGdoPrimalDfs(2,set_sz,CUT,node_mark,edge_mark,data);
  CHECKRVALG(rval,CLEANUP);

  /* test that the path is only within one orientation */
  cur_mark = node_mark[EGmengerGetEdgeData(path[0],os)->v];
  *orientation = cur_mark - 1;
  for( i = pathsz ; i-- ; )
  {
    pedge = EGmengerGetEdgeData(path[i],os);
    TESTG((rval=((node_mark[pedge->v]!=cur_mark) || 
                 (node_mark[GETOTHEREND(pedge)]!=cur_mark))), 
          CLEANUP, "Imposible!");
  }

  /* ending */
  CLEANUP:
  EGfreeList(CUT[0]);
  EGfreeList(CUT[1]);
  return rval;
}

/* ========================================================================= */
int EGgetDualCut( EGgreedyData_t*const data,
                  EGdualCut_t**const dualcut,
                  const int pset_sz,
                  const int*const pset,
                  unsigned char*const node_mark,
                  unsigned char*const edge_mark,
                  EGdGraphEdge_t**const dc_nodes)
{
  int rval = 0;
  register unsigned int i;
  unsigned cycle_sz = 0;
  EGlistNode_t *nit,*eit;
  EGdGraphNode_t*c_node;
  EGdGraphEdge_t*c_edge;
  const graphP G = data->G;
  EGdijkstraCost_t cutval = EGdijkstraToCost(0.0);
  graphNodeP pedge;
  /* reset marks */
  memset(dc_nodes,0,sizeof(EGdGraphEdge_t*)*data->norig_edges);
  memset(node_mark,0,sizeof(unsigned char)*data->norig_nodes);
  memset(edge_mark,0,sizeof(unsigned char)*data->norig_edges);
  /* mark all nodes in the cut */
  for( i = pset_sz ; i-- ; ) node_mark[pset[i]] = 1;
  /* loop through all edges in the bi-dual graph */
  for( nit = data->bdG->nodes->begin ; nit ;  nit = nit->next)
  {
    c_node = (EGdGraphNode_t*)nit->this;
    for(eit = c_node->out_edges->begin ; eit ; eit = eit->next)
    {
      c_edge = (EGdGraphEdge_t*)eit->this;
      pedge = EGmengerGetEdgeData(c_edge,os);
      if(edge_mark[pedge->id]) continue;
      if(node_mark[pedge->v] != node_mark[GETOTHEREND(pedge)])
      {
        edge_mark[pedge->id] = 1;
        cutval = EGdijkstraCostAdd( cutval, 
                                    EGdijkstraGetEdgeLength( c_edge, dij_os) );
        dc_nodes[cycle_sz++] = c_edge;
      }
    }
  }
  /* save the cut */
  *dualcut = EGnewDualCut(data->bdG->mem, cycle_sz);
  memcpy((*dualcut)->edges,dc_nodes,cycle_sz*sizeof(EGdGraphEdge_t*));
  (*dualcut)->value = cutval;
  return rval;
}

---