karger.c

Go to the documentation of this file.

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <time.h>
#include "karger.h"
#include "bc_zeit.h"
#include "eg_macros.h"

int max_numiter = 0;

#define MAXCKEYS 32768
#define MAXALPHA 4
#define MAXPSNODES 64
#define RFACTOR 1000

static int
   get_numiter(int nnode, int nedges, int alpha, double bound),
   create_cutset(int nnodes, int nps, int *psnodes, int *pssubset,
                 int nremain, int *edge, double *length, int *treeinf,
                 unsigned int *nkey, char *ckeys, int *cutset,
                 int *cutsize, double *cutweight, int *foundset),
  get_alphacut(int nnodes, int nedges, int *edge, double *length,
             int *treeinf, unsigned int *nkey, char *ckeys, int alpha,
             double bound, setlist **sets,
             int (*choose_edge)(int, int*, double*, double, void*, int*),
             void *process_info,
             int (*process_set)(int*, int, double, void*, setlist **)
             );

static void
    merge_components(int comp1, int comp2, int *treeinf, int *newnode);


/* nnodes, nedges, elist - graph with nnodes and nedges given in elist
   elen          - weight of edges
   alpha, bound  - want cuts with weight <= alpha * min-cut && <= bound
   maxtime, seed - maximum time, random seed (if 0, seed chosen from time)
   sets          - list of sets passed back
   choose_edge   - function to choose edge to be contracted, see choose_edge1
   process_info  - arbitrary info to be used in process_set
   process_set   - function to process and/or generated sets
*/
int karger(int nnodes, int nedges, int *elist, double *elen, int alpha,
           double bound, double maxtime, int seed, setlist **sets,
           int (*choose_edge)(int, int*, double*, double, void*, int*),
           void *process_info,
           int (*process_set)(int*, int, double, void*, setlist **)
           )
{
   int rval = PROCESS_OVER;
   int r = 0, i = 0;
   int *edge = NULL;
   double *length = NULL;
   int *treeinf = NULL;
   unsigned int *nkey = NULL, *nkeyp = NULL;
   char ckeys[MAXCKEYS];
   int numiter;
   double szeit = CCutil_zeit();

   if (nnodes <= 0 || nedges <= 0){
      fprintf(stderr,"Invalid # of nodes or # of edges\n");
      return PROCESS_ERROR;
   }
   if (elist == NULL || elen == NULL){
      fprintf(stderr,"Invalid edge or weight lists\n");
      return PROCESS_ERROR;
   }
   if (alpha < 0){
      fprintf(stderr,"Invalid alpha\n");
      return PROCESS_ERROR;
   }
   if (alpha > MAXALPHA){
      fprintf(stderr,"give alpha <= %d\n", MAXALPHA);
      return PROCESS_ERROR;
   }      
   if (bound <= 0){
      fprintf(stderr,"Invalid bound %f on cuts\n", bound);
      return PROCESS_ERROR;
   }
   if (sets == NULL){
      fprintf(stderr,"Null set list pointer\n");
      return PROCESS_ERROR;
   }
  if(seed == 0 ) seed = CCutil_zeit();
  srandom((unsigned int)seed);

   edge = MALLOC(2*nedges, int);
   length = MALLOC(nedges, double);
   treeinf = MALLOC(3*nnodes, int);
   nkey = MALLOC(nnodes, unsigned int);
   nkeyp = MALLOC(nnodes, unsigned int);
   if (edge == NULL || length == NULL || treeinf == NULL || nkey == NULL){
      fprintf(stderr,"Not enough mem for local arrays\n");
      rval = PROCESS_ERROR;
      goto CLEANUP;
   }

   /* copy the edges and lengths into local arrays */
   for (i=0; i<nedges; i++){
      edge[2*i] = elist[2*i];
      edge[2*i+1] = elist[2*i+1];
      length[i] = elen[i];
   }

   for (i=0; i<MAXCKEYS; i++)
      ckeys[i] = 0;
   for (i=0; i<nnodes; i++)
      nkeyp[i] = random();

   /*
   *sets = NULL;
   */
   numiter = get_numiter(nnodes, nedges, alpha, bound);
   for (r=0; r<numiter; r++){
      /* Initialize component info */
      for (i=0; i<nnodes; i++){
         treeinf[3*i] = i;
         treeinf[3*i+1] = 1;
         treeinf[3*i+2] = 0;

         /* assign random key to each node */
         nkey[i] = nkeyp[i];
      }

      if (maxtime < CCutil_zeit()-szeit){
         rval = PROCESS_OVER;
         goto CLEANUP;
      }

      rval = get_alphacut(nnodes, nedges, edge, length, treeinf, nkey, ckeys,
                          alpha, bound, sets, choose_edge, process_info,
                          process_set); 
      if (rval == PROCESS_OVER || rval == PROCESS_ERROR) 
      {
         goto CLEANUP;
      }
   }
   if(r == numiter) rval = 0;

CLEANUP:

   fprintf(stdout, "\nKarger: number of iterations: %d, time %lf\n", r, 
            CCutil_zeit()-szeit);

   IFFREE(edge, int);
   IFFREE(length, double);
   IFFREE(treeinf, int);
   IFFREE(nkey, unsigned int);
   IFFREE(nkeyp, unsigned int);

   return rval;

} /* end kargers */


int choose_edge1(int nremain, int *edge, double *length, double tweight,
                 void *process_info, int *e)
{
   int i, rval = 0;
   long maxwt;
   double twt, drandn;

   /* locate edge to be contracted */
   maxwt = (long)((tweight * (double)RFACTOR) + 0.5);
   drandn = (random() % maxwt) / RFACTOR;

   for (i=0, twt=0.0; i<nremain; i++){
      twt += length[i];
      if (drandn <= twt){
         *e = i;
         return rval;
      }
   }
   rval = 1;
   return rval;
}

int add_cut(int *cutset, int cutsize, double cutweight, void *process_info,
            setlist **sets)
{
   int i = 0;
   setlist *tset = NULL;
   int rval = PROCESS_INCOMPLETE;

   tset = (setlist *) calloc(1, sizeof(setlist));
   if (tset == NULL){
      fprintf(stderr, "No memory in add_cut\n");
      return PROCESS_ERROR;
   }
   tset->setv = (int *) calloc((size_t)cutsize, sizeof(int));
   if (tset->setv == NULL){
      fprintf(stderr, "No memory in process_tsp\n");
      IFFREE(tset, setlist);
      return PROCESS_ERROR;
   }
   tset->setn = cutsize;
   tset->cutval = cutweight;
   tset->next = *sets;
   for (i=0; i<cutsize; i++)
      tset->setv[i] = cutset[i];
   *sets = tset;

   return rval;
}

int karg_getprob (char *fname, int *nnodes, int *nedges, int **elist, 
                  double **elen)
{
    FILE *f = NULL;
    int i = 0;
    int n1 = 0, n2 = 0;
    double weight = 0;

    *elist = (int *) NULL;
    *elen  = (double *) NULL;

    if ((f = fopen (fname, "r")) == NULL) {
       fprintf (stderr, "Unable to open %s for input\n",fname);
       goto CLEANUP;
    }

    if (fscanf (f, "%d %d", nnodes, nedges) != 2) {
       fprintf (stderr, "File %s has invalid format\n",fname);
       goto CLEANUP;
    }
    printf ("Nodes: %d  Edges: %d\n", *nnodes, *nedges);
    fflush (stdout);

    *elist = MALLOC (2 * (*nedges), int);
    *elen  = MALLOC (*nedges, double);
    if (*elist == NULL || *elen == NULL) {
        fprintf (stderr, "out of memory in karg_getprob\n");
        goto CLEANUP;
    } 

    for (i = 0; i < *nedges; i++) {
       if (fscanf(f,"%d %d %lf", &n1, &n2, &weight) != 3) {
          fprintf (stderr, "%s has invalid input format\n",fname);
          goto CLEANUP;
       }
       if ((weight < 0) || (n1 == n2) || (n1 < 0) || (n2 < 0) ||
           (n1 >= *nnodes) || (n2 >= *nnodes)){
          fprintf (stderr, "%s has invalid input numbers\n",fname);
          goto CLEANUP;
       }
       (*elist)[2*i] = n1;
       (*elist)[2*i+1] = n2;
       (*elen)[i] = weight;
    }
    fclose (f);
    return 0;

CLEANUP:
    
    *nnodes = 0;
    *nedges = 0;
    IFFREE(*elist, int);
    IFFREE(*elen, double);
    if (f) fclose (f);
    return 1;
}

int process_cut(int *cutset, int cutsize, double cutweight, void *process_info,
                setlist **sets)
{
   int i = 0;
   setlist *tset = NULL;
   int rval = PROCESS_INCOMPLETE;
   cuts_info *cinf = (cuts_info *)process_info;

   /* If cinf.val or 0 is in set, then ignore */
   for (i=0; i<cutsize; i++)
      if (cutset[i] == 0 || cutset[i] == cinf->val)
         return rval;

   /* if cut has node 2 stop
   for (i=0; i<cutsize; i++)
      if (cutset[i] == 2)
         return PROCESS_OVER;
   */
   tset = (setlist *) calloc(1, sizeof(setlist));
   if (tset == NULL){
      fprintf(stderr, "No memory in add_cut\n");
      return PROCESS_ERROR;
   }
   tset->setv = (int *) calloc((size_t)cutsize, sizeof(int));
   if (tset->setv == NULL){
      fprintf(stderr, "No memory in process_tsp\n");
      IFFREE(tset, setlist);
      return PROCESS_ERROR;
   }
   tset->setn = cutsize;
   tset->cutval = cutweight;
   tset->next = *sets;
   for (i=0; i<cutsize; i++)
      tset->setv[i] = cutset[i];
   *sets = tset;

   return rval;
}

/***************************************/
/* From now on internal functions only */


/* modify this function to return whatever you want */
static int
get_numiter(int nnode, int nedges, int alpha, double bound)
{
   if (max_numiter == 0)
      return (alpha * nnode * (nnode - 1)) / 2;
   else
      return max_numiter;
}

static int
get_alphacut(int nnodes, int nedges, int *edge, double *length,
             int *treeinf, unsigned int *nkey, char *ckeys, int alpha,
             double bound, setlist **sets,
             int (*choose_edge)(int, int*, double*, double, void*, int*),
             void *process_info,
             int (*process_set)(int*, int, double, void*, setlist **)
             )
{
   int rval = PROCESS_INCOMPLETE;
   int ncontract = 0, nremain = 0, foundset = 0;
   int i = 0, j = 0, k = 0, e = 0;
   int newnode = 0, temp = 0;
   int comp1 = 0, comp2 = 0;
   int maxsets = 0, twoalpha = 0;
   int psnodes[MAXPSNODES],pssubset[MAXPSNODES], nps = 0;
   int *cutset = NULL, *setind = NULL, cutsize = 0;
   double temp1 = 0, tweight = 0;
   double cutweight = 0;

   if ((cutset = CALLOC(nnodes, int)) == NULL){
      fprintf(stderr, "No memory in get_alphacut()\n");
      return PROCESS_ERROR;
   }
   
   /* add up all edges */
   for (i=0, tweight=0.0; i<nedges; i++){
      tweight += length[i];
   }
   
   /* all edges remain, contracted 0 */
   nremain = nedges;
   ncontract = 0;

   /* need two X alpha nodes */
   twoalpha = alpha *2;
   
   /* contract till you have #nodes  = twoalpha left */
   for (ncontract = 0; ncontract < nnodes - twoalpha; ncontract++){
      if (nremain == 0)
         break;

      /* locate edge to be contracted */
      rval = choose_edge(nremain, edge, length, tweight, process_info, &e);
      if (rval) goto CLEANUP;

      if (e < 0 || e >= nremain){
         rval = PROCESS_ERROR;
         goto CLEANUP;
      }

      /* merge the two pseudo-nodes which form ends of the edges */
      comp1 = treeinf[ 3*edge[2*e] ];
      comp2 = treeinf[ 3*edge[2*e+1] ];
      merge_components(comp1, comp2, treeinf, &newnode);

      /* merge key */
      if (newnode == comp2)
         nkey[comp2] ^= nkey[comp1];
      else
         nkey[comp1] ^= nkey[comp2];

      /* move all edges lying within the pseudonode to the end */
      for (i=nremain-1; i>=0; i--){
         if (treeinf[ 3*edge[2*i] ] == treeinf[ 3*edge[2*i+1] ]){
            if (i != nremain-1){
               SWAP(edge[2*i], edge[2*(nremain-1)], temp);
               SWAP(edge[2*i+1], edge[2*(nremain-1)+1], temp);
               SWAP(length[i], length[nremain-1], temp1);
            }
            tweight -= length[nremain-1];
            nremain--;
         }
      }

   } /* end edges */

   if (nnodes - ncontract > MAXPSNODES){
      rval = PROCESS_ERROR;
      goto CLEANUP;
   }

   for (i=0, nps=0; i<nnodes; i++)
      if (treeinf[3*i+1] > 0) psnodes[nps++] = i;

   /* create indicator array */
   if ((setind = CALLOC(nnodes, int)) == NULL){
      fprintf(stderr, "No memory in get_alphacut()\n");
      rval = PROCESS_ERROR; goto CLEANUP;
   }

   /* If no edges remain, every node represents a cut */
   if (nremain == 0){
      for (i=0; i< nps; i++){
         pssubset[i] = 1;

         /* create cut-set */
         rval = create_cutset(nnodes, 1, &(psnodes[i]),&(pssubset[i]), nremain,
                              edge, length, treeinf, nkey, ckeys, cutset,
                              &cutsize, &cutweight, &foundset);
         if (rval) {rval = PROCESS_ERROR; goto CLEANUP; }
         if (foundset) continue;
         if (cutweight >= bound - TOLER) continue;

         /* process set */
         rval = process_set(cutset, cutsize, cutweight, process_info, sets);

         if (rval == PROCESS_ERROR || rval == PROCESS_OVER)
            goto CLEANUP;

         /* create negation */
         {
            //int k;

            for (k=0; k<nnodes; k++) setind[k] = 0;
            setind[psnodes[i]] = 1;
            for (k=0, cutsize = 0; k<nnodes; k++)
               if (!setind[treeinf[3*k]])
                  cutset[cutsize++] = k;
         }

         /* process negation */
         rval = process_set(cutset, cutsize, cutweight, process_info, sets);

         if (rval == PROCESS_ERROR || rval == PROCESS_OVER)
            goto CLEANUP;
      }
   }
   else{ /* nremain != 0 */
      /* generate sets using simple gray code */
      for (i=0; i<nps; i++)
         pssubset[i] = 0;

      maxsets = (1 << (nps-1)) - 1;
      for (i=0; i<maxsets; i++){
         for (j=0, k=1; j < nps; j++, k <<= 1){
            if ((k & i) == 0){
               pssubset[j] ^= 1;
               break;
            }
         }
         /* create cut-set */
         rval = create_cutset(nnodes, nps, psnodes, pssubset, nremain, edge,
                              length, treeinf, nkey, ckeys, cutset, &cutsize,
                              &cutweight, &foundset);
         if (rval) {rval = PROCESS_ERROR; goto CLEANUP; }
         if (foundset) continue;
         if (cutweight >= bound - TOLER) continue;

         /* process set */
         rval = process_set(cutset, cutsize, cutweight, process_info, sets);

         if (rval == PROCESS_ERROR || rval == PROCESS_OVER)
            goto CLEANUP;

         /* create negation */
         {
            //int k;

            for (k=0; k<nnodes; k++) setind[k] = 0;

            for (k=0; k<nps; k++)
               if (pssubset[k])
                  setind[psnodes[k]] = 1;
            for (k=0, cutsize = 0; k<nnodes; k++)
               if (!setind[treeinf[3*k]])
                  cutset[cutsize++] = k;
         }

         /* process negation */
         rval = process_set(cutset, cutsize, cutweight, process_info, sets);
         
         if (rval == PROCESS_ERROR || rval == PROCESS_OVER)
            goto CLEANUP;

      } /* end for all subsets */
   } /* nremain > 0 */

   rval = PROCESS_INCOMPLETE;

 CLEANUP:
   IFFREE(cutset, int);
   IFFREE(setind, int);
   return rval;

} /* End get cut */

static void
merge_components(int comp1, int comp2, int *treeinf, int *newnode)
{
   int cur, i, prev;
   int siz1, siz2, c1, c2;

   /* merge smaller comp */
   siz1 = treeinf[3*comp1+1];
   siz2 = treeinf[3*comp2+1];

   if (siz1 <= siz2){
      c1 = comp1;
      c2 = comp2;
   }
   else{
      c2 = comp1;
      c1 = comp2;
      siz1 = siz2;
   }

   for (i=0,cur=c1; i<siz1; i++){
      treeinf[3*cur] = c2;
      prev = cur;
      cur = treeinf[3*cur+2];
   }
   treeinf[3*c2+1] += siz1;
   treeinf[3*c1+1] = 0;
   treeinf[3*prev+2] = treeinf[3*c2+2];
   treeinf[3*c2+2] = c1;
   *newnode = c2;
}


/* This assumes that the function will be called for only one side of a cut,
   not for it's negation.
*/
static int
create_cutset(int nnodes, int nps, int *psnodes, int *pssubset,
              int nremain, int *edge, double *length, int *treeinf,
              unsigned int *nkey, char *ckeys, int *cutset,
              int *cutsize, double *cutweight, int *foundset)
{
   int i = 0,tsiz = 0;
   unsigned int rkey = 0;
   int *setind = NULL, csiz = 0, stsiz = 0;
   int rval = 0;

   setind = CALLOC(nnodes, int);
   if (setind == NULL){
      fprintf(stderr, "No memory in create_cutset()\n");
      rval = 1;
      goto CLEANUP;
   }

   for (i=0, stsiz=0, csiz = 0; i<nps; i++){
      if (pssubset[i]){
         csiz += treeinf[3*psnodes[i]+1];
         setind[psnodes[i]] = 1;
         stsiz++;
      }
   }

   if (nps < 1 || stsiz < 1 || csiz < 1){
      fprintf(stderr, "Invalid input to process_set()\n");
      rval = 1;
      goto CLEANUP;
   }

   /* form the key for the smaller side of the cut and store */
   if (nnodes - csiz > csiz){
      for (i=0, rkey=0; i<nps; i++)
         if (pssubset[i])
            rkey ^= nkey[psnodes[i]];
   }
   else{
      for (i=0, rkey=0; i<nps; i++)
         if (!pssubset[i])
            rkey ^= nkey[psnodes[i]];
   }
   if (ckeys[rkey % MAXCKEYS] == 0){
      *foundset = 0;
      ckeys[rkey % MAXCKEYS] = 1;
   }
   else{ /* set found with high prob */
      *foundset = 1;
      goto CLEANUP;
   }

   for (i=0, (*cutweight)=0.0; i<nremain; i++)
      if (setind[ treeinf[ 3*edge[2*i] ]] != setind[treeinf[ 3*edge[2*i+1]]])
         (*cutweight) += length[i];

   *cutsize = csiz;
   for (i=0, csiz=0; i<nnodes; i++)
      if (setind[treeinf[3*i]])
         cutset[csiz++] = i;

   if (csiz != *cutsize){
      fprintf(stderr, "Some mistake in input to create_cutset\n");
      rval = 1;
   }

 CLEANUP:
   IFFREE(setind, int);
   return rval;
}

---