graphdual.c

Go to the documentation of this file.

#include "graphdual.h"
#ifndef GRD_TEST
#define GRD_TEST 1
#endif

/* this is an intermediate structure to build our graph */
typedef struct
{
  double value;
  int tail;
  int head;
  int index;
}
localEdge_t;

typedef const localEdge_t *clep_t;
typedef const clep_t *cclep_t;

static int localEdgeCompare (const void *p1,
                             const void *p2)
{
  return ((*((cclep_t) (p1)))->value < (*((cclep_t) (p2)))->value) ? 1 : -1;
}

/* to sort edge indices according to weight (indicated by x). 
   assumes that the indices array has already been allocated
   and is of size at least nedges */
static int sortEdgeIndices (int nedges,
                            int *edges,
                            int *indices,
                            double *x)
{

  int n;
  localEdge_t **edgeArray;

  /* now we need to create the edge array */
  edgeArray = (localEdge_t **) malloc (sizeof (localEdge_t *) * nedges);
  for (n = 0; n < nedges; n++)
  {
    edgeArray[n] = (localEdge_t *) malloc (sizeof (localEdge_t));
    edgeArray[n]->value = x[n];
    edgeArray[n]->head = edges[2 * n];
    edgeArray[n]->tail = edges[2 * n + 1];
    edgeArray[n]->index = indices[n];
  }

  /* now we sort the local edges */
  if (nedges > 1)
    qsort (edgeArray, (unsigned) nedges, sizeof (localEdge_t *),
           localEdgeCompare);
  EXIT (edgeArray[0]->value < edgeArray[nedges - 1]->value,
        "edgeArray is not well sorted");

  for (n = 0; n < nedges; n++)
  {
    indices[n] = edgeArray[n]->index;
    free (edgeArray[n]);
  }

  free (edgeArray);

  return 0;

}

/* to sort edges according to weight (indicated by x) */
static int sortEdges (int nedges,
                      int *edges,
                      double *x)
{

  int n;
  localEdge_t **edgeArray;

  /* now we need to create the edge array */
  edgeArray = (localEdge_t **) malloc (sizeof (localEdge_t *) * nedges);
  for (n = 0; n < nedges; n++)
  {
    edgeArray[n] = (localEdge_t *) malloc (sizeof (localEdge_t));
    edgeArray[n]->value = x[n];
    edgeArray[n]->head = edges[2 * n];
    edgeArray[n]->tail = edges[2 * n + 1];
    edgeArray[n]->index = n;
  }

  /* now we sort the local edges */
  if (nedges > 1)
    qsort (edgeArray, (unsigned) nedges, sizeof (localEdge_t *),
           localEdgeCompare);
  EXIT (edgeArray[0]->value < edgeArray[nedges - 1]->value,
        "edgeArray is not well sorted");

  for (n = 0; n < nedges; n++)
  {
    edges[2 * n] = edgeArray[n]->head;
    edges[2 * n + 1] = edgeArray[n]->tail;
    x[n] = edgeArray[n]->value;
    free (edgeArray[n]);
  }

  free (edgeArray);

  return 0;

}

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

/* generate the dual graph and dual planar embedding, and dual planar embeding
 * with links to the actual outgoing edges in the bidirected dual. */
EGdGraph_t *getDualBoyer (graphP G,
                          int nnodes,
                          int nedges,
                          double *weigh,
                          EGlist_t *** dembed,
                          EGmemPool_t * localmem,
                          EGlist_t *** lembed)
{

  /* local variables */
  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)
  };
  EGdGraph_t *bdG = 0;
  EGdGraphEdge_t *e;
  EGdGraphNode_t **node_array;
  EGmengerNodeData_t *ndata;
  EGmengerEdgeData_t *edata;
  int nroot,
    nnexti,
    i,
    rval,
    curface,
    nnextii;
  graphNodeP eproot,
    eroot,
    epnexti,
    epnextii;
  unsigned short *e_old,
   *e_mark;
  const int nFaces = nedges - nnodes + 2;
  unsigned char *face_mark;

  /* basic set-up */
  ADVTESTL (GD_LEVEL, nnodes != G->N, 0, "number nodes in graph and "
            "indicated doesn't match %d != %d", nnodes, G->N);
  ADVTESTL (GD_LEVEL, nedges != G->M, 0, "number edges in graph and "
            "indicated doesn't match %d != %d", nedges, G->M);
  e_old = EGmemPoolSMalloc (localmem, unsigned short,
                            G->M);
  e_mark = EGmemPoolSMalloc (localmem, unsigned short,
                             G->M);
  face_mark = EGsMalloc (unsigned char,
                         nFaces);
  memset (face_mark, 0, sizeof (unsigned char) * nFaces);
  node_array = EGmemPoolSMalloc (localmem, EGdGraphNode_t *, nFaces);
  *dembed = EGmemPoolSMalloc (localmem, EGlist_t *, nFaces);
  if(lembed) *lembed = EGmemPoolSMalloc (localmem, EGlist_t *, nFaces);
  bdG = EGnewDGraph (localmem);

  i = gp_Embed (G, EMBEDFLAGS_PLANAR);
  rval = gp_SortVertices (G);
  EXITRVAL (rval);
  if (i != OK)
    return 0;

  /* now we have the embedding, we only need to count how 
   * many faces there are there and to create bdG, we start
   * setting mark and old to 0 for both nodes and edges */
  for (nroot = 0; nroot < G->N; nroot++)
  {
    if (ISNODEEMPTY (nroot, G))
      continue;
    eproot = GETROOTEDGE (nroot, G);
    e_mark[eproot->id] = 0;
    e_old[eproot->id] = 0;
    /* we set all mark edges in the adjacency list to zero */
    for (epnexti = GETNEXTEDGE (eproot, nroot, G); epnexti != eproot;
         epnexti = GETNEXTEDGE (epnexti, nroot, G))
    {
      e_mark[epnexti->id] = 0;
      e_old[epnexti->id] = 0;
    }
  }

  /* ======================================================================= */
  /* FIRST STEP: SET THE MARCS IN OLD AND IN MARK */
  /* now we look at every edge pair, but discard if we have 
   * already traverse the edge in that direction, if we have 
   * traversed the edge in the direction end0->end1, then 
   * mark = nFace, if we have traversed the edge in the direction 
   * end1->end0, then old = nFace */
  curface = 1;
  /* loop through all nodes */
  for (nroot = 0; nroot < G->N; nroot++)
  {
    /* eproot point to the first edge in the embbedding list */
    eproot = GETROOTEDGE (nroot, G);
    epnexti = eproot;
    do
    {
      /* loop through neighbours */
      if (((unsigned) (epnexti - G->G)) & 1U)
      {
        /* if we already traversed this edge in this direction we are done */
        if (e_mark[epnexti->id])
        {
          epnexti = GETNEXTEDGE (epnexti, nroot, G);
          continue;
        }
      }
      else
      {
        /* if we already traversed this edge in this direction we are done */
        if (e_old[epnexti->id])
        {
          epnexti = GETNEXTEDGE (epnexti, nroot, G);
          continue;
        }
      }
      /* we select nnextii as the next end point of this edge */
      nnextii = epnexti->v;

      /* now we loop through the edges of this face */
      epnextii = epnexti;
      nnexti = nroot;
      eroot = epnextii;
      do
      {                         /* loop through the face */
        if (((unsigned) (epnextii - G->G)) & 1U)
        {
          e_mark[epnextii->id] = curface;
          i = e_old[epnextii->id];
        }
        else
        {
          e_old[epnextii->id] = curface;
          i = e_mark[epnextii->id];
        }
        /* we update epnextii, eroot, nnexti and nnextii */
        epnextii = GETROOTEDGE (nnextii, G);
        /* find the epnextii that points to eroot */
        while (epnextii->v != nnexti)
          epnextii = GETNEXTEDGE (epnextii, nnextii, G);
        /* move according to the embedding of edge */
        epnextii = GETNEXTEDGE (epnextii, nnextii, G);
        eroot = epnextii;
        nnexti = nnextii;
        nnextii = eroot->v;
      } while ((epnextii != epnexti) || (nroot != nnexti));
      /* end loop through face */
      /* at this stage we have finish a face, and now we create a new one */
      curface++;
      EXIT (curface - 1 > nFaces, "Number of faces is wrong nFaces %u "
            "curface %u", nFaces, curface - 1);
      /*update epnexti */
      epnexti = GETNEXTEDGE (epnexti, nroot, G);
    } while (epnexti != eproot);  /* end loop through neighbours */
  }                             /* end loop through all nodes */

  /* ======================================================================= */
  /* SECOND PHASE, CREATE DUAL GRAPH AND DUAL EMBEDING */
  /* now we look at every edge pair, but discard if we have 
   * already traverse the edge in that direction, if we have 
   * traversed the edge in the direction end0->end1, then 
   * mark = nFace, if we have traversed the edge in the direction 
   * end1->end0, then old = nFace */
  /* create the nodes and embed lists */
  for (i = 0; i < nFaces; i++)
  {
    ndata = EGnewMengerNodeData (localmem);
    node_array[i] = EGdGraphAddNode (bdG, ndata);
    (*dembed)[i] = EGnewList (localmem);
    if(lembed) (*lembed)[i] = EGnewList (localmem);
  }
  /* loop through all nodes */
  for (nroot = 0; nroot < G->N; nroot++)
  {
    /* eproot point to the first edge in the embbedding list */
    eproot = GETROOTEDGE (nroot, G);
    epnexti = eproot;
    do
    {
      /* loop through neighbours */
      if (((unsigned) (epnexti - G->G)) & 1U)
      {
        /* if we already traversed this edge in this direction we are done */
        curface = e_mark[epnexti->id];
        if (face_mark[curface - 1])
        {
          epnexti = GETNEXTEDGE (epnexti, nroot, G);
          continue;
        }
      }
      else
      {
        /* if we already traversed this edge in this direction we are done */
        curface = e_old[epnexti->id];
        if (face_mark[curface - 1])
        {
          epnexti = GETNEXTEDGE (epnexti, nroot, G);
          continue;
        }
      }
      face_mark[curface - 1] = 1;
      /* we select nnextii as the next end point of this edge */
      nnextii = epnexti->v;

      /* now we loop through the edges of this face */
      epnextii = epnexti;
      nnexti = nroot;
      eroot = epnextii;
      do
      {                         /* loop through the face */
        if (((unsigned) (epnextii - G->G)) & 1U)
        {
          i = e_old[epnextii->id];
        }
        else
        {
          i = e_mark[epnextii->id];
        }
        /* now we add the corresponding edge to the graph and the embed list */
        EGlistPushBack ((*dembed)[curface - 1], (void *) (epnextii->id));
        edata = EGnewMengerEdgeData (localmem);
        e = EGdGraphAddEdge (bdG, edata,
                             node_array[i - 1], node_array[curface - 1]);
        EGmengerSetEdgeCost (e, os, EGdijkstraToCost (weigh[epnextii->id]));
        EGmengerSetEdgeData (e, os, epnextii);
        EGmengerSetEdgeIsInSolution (e, os, 0);
        if(lembed) EGlistPushBack((*lembed)[curface - 1], e);
        /* we update epnextii, eroot, nnexti and nnextii */
        epnextii = GETROOTEDGE (nnextii, G);
        /* find the epnextii that points to eroot */
        while (epnextii->v != nnexti)
          epnextii = GETNEXTEDGE (epnextii, nnextii, G);
        /* move according to the embedding of edge */
        epnextii = GETNEXTEDGE (epnextii, nnextii, G);
        eroot = epnextii;
        nnexti = nnextii;
        nnextii = eroot->v;
      } while ((epnextii != epnexti) || (nroot != nnexti)); /* end loop through face */
      /*update epnexti */
      epnexti = GETNEXTEDGE (epnexti, nroot, G);
    } while (epnexti != eproot);  /* end loop through neighbours */
  }                             /* end loop through all nodes */

  /* ending */
  EGmemPoolSFree (node_array, EGdGraphNode_t *, nFaces, localmem);
  EGmemPoolSFree (e_old, unsigned short,
                  G->M,
                  localmem);
  EGmemPoolSFree (e_mark, unsigned short,
                  G->M,
                  localmem);
  EGfree (face_mark);
  return bdG;
}

/* generate a planar subgraph from a graph by building a maximum weight spanning
 * tree and then eliminating edges that make it non planar */
int DPedgeEliminationHeuristic (int nnodes,
                                int nedges,
                                int *edges,
                                double *weight,
                                int *nplanar_edges,
                                int *planar_edges,
                                double *planar_weight,
                                int *nelim_indices,
                                int *elim_indices)
{

#if DISPLAY_MODE
  fprintf (stdout, "EEH: Initiating Edge Elimination Heuristic.\n");
  fflush (stdout);
#endif

  int rval;

  rval = sortEdges (nedges, edges, weight);
  CHECKRVAL (rval);

#if DISPLAY_MODE
  fprintf (stdout, "EEH: Executing binary search.\n");
  fflush (stdout);
#endif

  rval =
    DPbinPlanarizeBoyer (nnodes, nedges, 0, edges, weight, nplanar_edges,
                         planar_edges, planar_weight, elim_indices);
  CHECKRVAL (rval);

  *nelim_indices = nedges - *nplanar_edges;

#if DISPLAY_MODE
  fprintf (stdout, "EEH: Finished.\n");
  fflush (stdout);
#endif

  return 0;

}

int isPlanarBoyer (int nnodes,
                   int nedges,
                   int *edges)
{

  int rval,
    is_G_planar_b;
  graphP G;

  G = gp_New ();
  rval = cook2boyerGraph (G, nnodes, nedges, edges);
  EXITRVAL (rval);
  is_G_planar_b = gp_Embed (G, EMBEDFLAGS_PLANAR);
  gp_Free (&G);

  if (is_G_planar_b == OK)
    return 1;

  return 0;

}

int isPlanarOrMinorBoyer (int nnodes,
                          int nedges,
                          int *edges,
                          int *nmedges,
                          int *medges)
{

  int rval,
    is_G_planar_b,
    nmnodes;
  graphP G;

  G = gp_New ();
  rval = cook2boyerGraph (G, nnodes, nedges, edges);
  EXITRVAL (rval);
  is_G_planar_b = gp_Embed (G, EMBEDFLAGS_PLANAR);

  if (is_G_planar_b == OK)
  {
    gp_Free (&G);
    return 1;
  }

  gp_SortVertices (G);
  extractCookEdgeIds (G, &nmnodes, nmedges, medges);
  gp_Free (&G);

  return 0;

}

/* This function works with two graphs G1 and G2 spanning a common set of
   vertices. It assumes that G2 is a planar graph, and that its edge set 
   is strictly contained in that of G1. It finds a graph G3 in between
   G1 and G2 that is planar by means of binary search and edge additions.
  
   nnodes: number of nodes in G1 and G2.
   nedges1: number of edges of G1.
   nedges2: number of edges of G2.
   edges: array of edges in G1. It is assumed that the first nedges2 edges
          of the array are the edges of G2. In Cook format.
   nedges3: number of edges in G3 (output).
   edges3: array with edges in G3 in Cook format (output).
   elim_indices: array with the indices of those edges eliminated from edges. */
int DPbinPlanarizeBoyer (int nnodes,
                         int nedges1,
                         int nedges2,
                         int *edges,
                         double *weigh,
                         int *nedges3,
                         int *edges3,
                         double *weigh3,
                         int *elim_indices)
{

  int nelim = 0,
    i,
    rval,
    current_edge,
    is_G_planar_b,
    top,
    bottom,
    middle,
    ltop,
    lbottom,
    lmiddle;
  graphP G;

#if DISPLAY_MODE
  fprintf (stdout, "\rBIN: Binary Search Heuristic.\n");
  fflush (stdout);
#endif

  /* Make sure that the graph is not-planar before starting                   */
  G = gp_New ();
  rval = cook2boyerGraph (G, nnodes, nedges1, edges);
  EXITRVAL (rval);
  is_G_planar_b = gp_Embed (G, EMBEDFLAGS_PLANAR);
  gp_Free (&G);

  /* If the graph is planar, return                                           */
  if (is_G_planar_b == OK)
  {
#if DISPLAY_MODE
    fprintf (stdout, "\rBIN: Graph is planar.\n");
    fflush (stdout);
#endif
    *nedges3 = 0;
    return 0;
  }

  /* G3 must start as a graph identical to G2                                 */

  *nedges3 = nedges2;

  for (i = 0; i < nedges2; i++)
  {
    edges3[2 * i] = edges[2 * i];
    edges3[2 * i + 1] = edges[2 * i + 1];
    weigh3[i] = weigh[i];
  }
  top = nedges2;
  bottom = nedges1;
  middle = top + (bottom - top) / 2;

  while ((bottom - top) > 0)
  {

    ltop = top;
    lbottom = bottom;
    lmiddle = ltop + (lbottom - ltop) / 2;

    while (lmiddle != ltop)
    {

      /* check if adding edges ltop to lmiddle makes graph planar             */

      current_edge = (*nedges3);
      for (i = ltop; i < lmiddle; i++)
      {
        edges3[2 * (*nedges3)] = edges[2 * i];
        edges3[2 * (*nedges3) + 1] = edges[2 * i + 1];
        weigh3[*nedges3] = weigh[i];
        (*nedges3)++;
      }

      /* check if adding it makes the graph non-planar                        */
      G = gp_New ();
      rval = cook2boyerGraph (G, nnodes, (*nedges3), edges3);
      EXITRVAL (rval);
      is_G_planar_b = gp_Embed (G, EMBEDFLAGS_PLANAR);
      gp_Free (&G);

      /* If the graph is planar, keep the edges and move the ltop down.       */
      if (is_G_planar_b == OK)
      {
        ltop = lmiddle;
        lmiddle = ltop + (lbottom - ltop) / 2;
        continue;
      }

      /* If the graph is not planar, remove the edges and move lbottom up.    */
      else
      {
        (*nedges3) = current_edge;
        lbottom = lmiddle;
        lmiddle = ltop + (lbottom - ltop) / 2;
        continue;
      }

    }

    /* We found the edge that needs to be eliminated */
    top = lbottom;
    elim_indices[nelim++] = lmiddle;

#if DISPLAY_MODE
    fprintf (stdout, "\rBIN: Removing edge (%6d %6d) [%lf]\n",
             edges[2 * lmiddle], edges[2 * lmiddle + 1], weigh[lmiddle]);
    fflush (stdout);
#endif

    /* If we add the rest of the edges, is the graph planar? */

    current_edge = (*nedges3);
    for (i = top; i < bottom; i++)
    {
      edges3[2 * (*nedges3)] = edges[2 * i];
      edges3[2 * (*nedges3) + 1] = edges[2 * i + 1];
      weigh3[(*nedges3)] = weigh[i];
      (*nedges3)++;
    }

    /* check if adding it makes the graph non-planar                        */
    G = gp_New ();
    rval = cook2boyerGraph (G, nnodes, (*nedges3), edges3);
    EXITRVAL (rval);
    is_G_planar_b = gp_Embed (G, EMBEDFLAGS_PLANAR);
    gp_Free (&G);

    /* If so, break. We are done. */
    if (is_G_planar_b == OK)
      break;

    /* Else, put the edges back and search for the next bad edge */

    else
      (*nedges3) = current_edge;

  }

#if DISPLAY_MODE
  fprintf (stdout, "\rBIN: Eliminated %d edges.\n", nelim);
  fprintf (stdout, "BIN: Completed planar sub-graph search.\n");
  fflush (stdout);
#endif

  return 0;

}

/* for very-safe shrinking... only forks */

int DPgetTrivialNodesToContract (int nnodes,
                                 int nedges,
                                 int *edges,
                                 double *weight,
                                 int *node_1,
                                 int *node_2)
{

  int i;
  int *degree;

  *node_1 = -1;
  *node_2 = -1;

  degree = (int *) calloc ((size_t) nnodes, sizeof (int));

  for (i = 0; i < nedges; i++)
    if (weight[i] > 0.99)
    {
      degree[edges[2 * i]] += 1;
      degree[edges[2 * i + 1]] += 1;
    }

  for (i = 0; i < nedges; i++)
    if (degree[edges[2 * i]] == 2 || degree[edges[2 * i + 1]] == 2)
    {
      *node_1 = edges[2 * i];
      *node_2 = edges[2 * i + 1];
      break;
    }

  //fprintf(stdout, "edge %d %d\n", *node_1, *node_2);

  free (degree);

  return 0;

}

int DPgetNonMinorNodesToContract (int nnodes,
                                  int nedges,
                                  int *edges,
                                  double *weight,
                                  int *node_1,
                                  int *node_2)
{

  *node_1 = -1;
  *node_2 = -1;

  int rval;
  int i,
    nmedges = 0,
   *medges,
   *degree,
    nmnodes = 0,
    max_degree = -1;

  medges = (int *) malloc (sizeof (int) * nedges);
  degree = (int *) malloc (sizeof (int) * nnodes);

  memset (degree, 0, sizeof (int) * nnodes);

  /*
   * {
   * rval = isPlanarOrMinorBoyer(nnodes, nedges, edges, &nmedges, medges);
   * 
   * if (rval)
   * goto CLEANUP;
   * }
   */

  {
    rval = isPlanarBoyer (nnodes, nedges, edges);
    if (rval)
      goto CLEANUP;

    rval = sortEdges (nedges, edges, weight);
    if (rval)
      goto CLEANUP;

    rval = DPgetBinMinor (nnodes, nedges, edges, &nmedges, medges, nedges);
    if (rval)
      goto CLEANUP;
  }

  max_degree = 0;

  for (i = 0; i < nmedges; i++)
  {
    degree[edges[2 * medges[i]]] += 1;
    degree[edges[2 * medges[i] + 1]] += 1;

    if (degree[edges[2 * medges[i]]] == 1)
      nmnodes++;
    if (degree[edges[2 * medges[i]]] > max_degree)
      max_degree = degree[edges[2 * medges[i]]];
    if (degree[edges[1 + 2 * medges[i]]] == 1)
      nmnodes++;
    if (degree[edges[1 + 2 * medges[i]]] > max_degree)
      max_degree = degree[edges[1 + 2 * medges[i]]];
  }

  /*
   * {
   * 
   * int n0 = 0, n2=0, n3=0, n4=0, nX=0;
   * 
   * for(i=0; i < nnodes; i++)
   * {
   * if ( degree[i] == 0 ) n0++;  
   * else if ( degree[i] == 2 ) n2++;
   * else if (degree[i] == 3) n3++;
   * else if (degree[i] == 4) n4++;
   * else nX++;
   * }
   * 
   * fprintf(stdout, "n0 = %d, n2 = %d, n3 = %d, n4 = %d, nX = %d\n", n0, n2, n3, n4, nX);
   * 
   * }
   */


#if 0 == 1
  /* first, try to contract 1-paths */
  for (i = 0; i < nmedges; i++)
    if (degree[edges[2 * medges[i]]] == 2
        && degree[edges[2 * medges[i] + 1]] == 2)
      if (weight[medges[i]] > 0.999)
      {
        *node_1 = edges[2 * medges[i]];
        *node_2 = edges[2 * medges[i] + 1];
        fprintf (stdout, "contracted strong path edge of weight 1\n");
        goto DISPLAY;
      }

  /* next, try to contract edges of weight = 1 */
  for (i = 0; i < nmedges; i++)
    if (degree[edges[2 * medges[i]]] == 2
        || degree[edges[2 * medges[i] + 1]] == 2)
      if (weight[medges[i]] > 0.999)
      {
        *node_1 = edges[2 * medges[i]];
        *node_2 = edges[2 * medges[i] + 1];
        fprintf (stdout, "contracted path edge of weight 1\n");
        goto DISPLAY;
      }
#endif

  /* find one edges that are unrelated to the minor */
  for (i = 0; i < nedges; i++)
    if (weight[i] > 0.99)
      if (!degree[edges[2 * i]] || !degree[edges[2 * i + 1]])
      {
        *node_1 = edges[2 * i];
        *node_2 = edges[2 * i + 1];
        fprintf (stdout, "contracted unrelated edge of weight 1\n");
        goto DISPLAY;
      }

#if 0 == 1
  /* finally, contract minor-paths */
  for (i = 0; i < nmedges; i++)
    if (degree[edges[2 * medges[i]]] == 2
        || degree[edges[2 * medges[i] + 1]] == 2)
    {
      *node_1 = edges[2 * medges[i]];
      *node_2 = edges[2 * medges[i] + 1];
      fprintf (stdout, "contracted path edge\n");
      goto DISPLAY;
    }
#endif

DISPLAY:

  ADVTESTL (0, max_degree != 3
            && max_degree != 4, 1, "Error! unidentified minor");

  if (*node_1 == -1)
    saveSubGraph ("contracted_to_minor.x", nnodes, nmedges, edges, medges,
                  weight);

CLEANUP:

#if DISPLAY_MODE
  if (max_degree == 3)
    fprintf (stdout,
             "MINOR_CONTRACT: Found K3,3 minor (m=%d, n=%d). Nodes to contract: %d and %d.\n",
             nmedges, nmnodes, *node_1, *node_2);
  else if (max_degree == 4)
    fprintf (stdout,
             "MINOR_CONTRACT: Found   K5 minor (m=%d, n=%d). Nodes to contract: %d and %d.\n",
             nmedges, nmnodes, *node_1, *node_2);
  else if (max_degree == -1)
    fprintf (stdout,
             "MINOR_CONTRACT: Graph is planar. No nodes to contract.\n");
  else
    fprintf (stdout, "MINOR_CONTRACT: Error! Error!\n");
  fflush (stdout);
#endif

  free (medges);
  free (degree);

  return 0;

}

int DPgetMinorNodesToContract (int nnodes,
                               int nedges,
                               int *edges,
                               int *node_1,
                               int *node_2)
{

  *node_1 = -1;
  *node_2 = -1;

  int rval,
    nfound = 0;
  int i,
    nmedges = 0,
   *medges,
   *degree,
    nmnodes = 0;

  medges = (int *) malloc (sizeof (int) * nedges);
  degree = (int *) malloc (sizeof (int) * nnodes);

  memset (degree, 0, sizeof (int) * nnodes);

  rval = isPlanarOrMinorBoyer (nnodes, nedges, edges, &nmedges, medges);

  if (rval)
    goto CLEANUP;

  for (i = 0; i < nmedges; i++)
  {
    degree[edges[2 * medges[i]]] += 1;
    degree[edges[2 * medges[i] + 1]] += 1;
    if (degree[edges[2 * medges[i]]] == 1)
      nmnodes++;
    if (degree[edges[1 + 2 * medges[i]]] == 1)
      nmnodes++;
  }

  for (i = 0; i < nnodes; i++)
  {
    if (degree[i] > 2 && !nfound)
    {
      *node_1 = i;
      nfound = 1;
    }
    else if (degree[i] > 2)
    {
      *node_2 = i;
      break;
    }
  }

  {

    int error = 0;

    if (degree[*node_1] == 3 && degree[*node_2] != 3)
      error = 1;
    else if (degree[*node_1] == 4 && degree[*node_2] != 4)
      error = 1;
    else if (degree[*node_1] != 3 && degree[*node_1] != 4)
      error = 1;

    if (error)
      for (i = 0; i < nnodes; i++)
        if (degree[i] != 0 && degree[i] != 2)
          fprintf (stderr, "degree[%d] = %d\n", i, degree[i]);
  }

  ADVTESTL (0, *node_1 == -1
            || *node_2 == -1, 1, "Failed to find nodes to contract");
  ADVTESTL (2, degree[*node_1] == 3
            && degree[*node_2] != 3, 1, "Error! unidentified minor");
  ADVTESTL (2, degree[*node_1] == 4
            && degree[*node_2] != 4, 1, "Error! unidentified minor");
  ADVTESTL (2, degree[*node_1] != 3
            && degree[*node_1] != 4, 1, "Error! unidentified minor");

#if DISPLAY_MODE
  if (degree[*node_1] == 3)
    fprintf (stdout,
             "MINOR_CONTRACT: Found K3,3 minor (n=%d, m=%d). Nodes to contract: %d and %d.\n",
             nmedges, nmnodes, *node_1, *node_2);
  else if (degree[*node_1] == 4)
    fprintf (stdout,
             "MINOR_CONTRACT: Found   K5 minor (n=%d, m=%d). Nodes to contract: %d and %d.\n",
             nmedges, nmnodes, *node_1, *node_2);
#endif

CLEANUP:

  free (medges);
  free (degree);

  return 0;

}

int DPfindBadEdgeK (int nnodes,
                    int nedges,
                    int *edges,
                    double *weight,
                    int k)
{

  int rval,
    bad_edge = -1;
  int who,
    i,
    nmedges = 0,
   *medges,
   *kedges = 0,
   *kindices = 0;
  double *kweight = 0,
    rvalue;

  if (k == 1)
    return (DPfindBadEdge (nnodes, nedges, edges, weight));
  else if (k < 1)
    return -1;

  medges = (int *) malloc (sizeof (int) * nedges);

  rval = isPlanarOrMinorBoyer (nnodes, nedges, edges, &nmedges, medges);

  //fprintf(stdout, "DPfindBadEdgeK: nmedges = %d\n", nmedges);

  if (rval)
    goto CLEANUP;

  if (k > nmedges)
    k = nmedges;

  kedges = (int *) malloc (sizeof (int) * 2 * nmedges);
  kweight = (double *) malloc (sizeof (double) * nmedges);
  kindices = (int *) malloc (sizeof (int) * nmedges);

  for (i = 0; i < nmedges; i++)
  {
    kedges[2 * i] = edges[2 * (medges[i])];
    kedges[2 * i + 1] = edges[2 * (medges[i]) + 1];
    kweight[i] = -1 * weight[medges[i]];
    kindices[i] = i;
  }

  sortEdgeIndices (nmedges, kedges, kindices, kweight);

  rvalue = random ();
  rvalue /= EGRAND_MAX;

  who = (int) (k * rvalue);

  if (who == k)
    who = k - 1;

  bad_edge = medges[kindices[who]];

CLEANUP:

  free (medges);

  if (kedges)
    free (kedges);
  if (kweight)
    free (kweight);
  if (kindices)
    free (kindices);

  return bad_edge;

}

int DPfindBadEdge (int nnodes,
                   int nedges,
                   int *edges,
                   double *weight)
{

  int rval;
  int i,
    nmedges = 0,
   *medges,
    min_ind = -1;
  double min = 200000000.0;

  medges = (int *) malloc (sizeof (int) * nedges);

  rval = isPlanarOrMinorBoyer (nnodes, nedges, edges, &nmedges, medges);

  if (rval)
    goto CLEANUP;

  for (i = 0; i < nmedges; i++)
    if (weight[medges[i]] < min)
    {
      min = weight[medges[i]];
      min_ind = medges[i];
    }

CLEANUP:

  free (medges);

  return min_ind;

}

/* This function works with a graph G1. It finds a graph G2 contained in G1
   which is planar by eliminating edges in Kuratowski minors.
  
   nnodes: number of nodes in G1.
   nedges1: number of edges of G1.
   edges: array of edges in G1. (Cook format).
   nedges3: number of edges in G3 (output).
   edges3: array with edges in G3 in Cook format (output).
   elim_indices: array with the indices of those edges eliminated from edges. */
int DPfastEdgeEliminationHeuristic (int nnodes,
                                    int nedges1,
                                    int *edges,
                                    double *weight,
                                    int *nedges2,
                                    int *edges2,
                                    double *weight2,
                                    int *nelim,
                                    int *elim_indices,
                                    int k_param)
{

  int i,
    iswap,
    bad_edge;
  //int bad_edge_2;
  int *indices = 0;

  double dswap;

#if DISPLAY_MODE
  fprintf (stdout, "\rSEE: Simple Edge Elimination Heuristic.\n");
  fflush (stdout);
#endif

  *nelim = 0;

  indices = (int *) malloc (sizeof (int) * nedges1);

  /* Copy G1 unto G2                                                          */

  *nedges2 = nedges1;
  for (i = 0; i < nedges1; i++)
  {
    edges2[2 * i] = edges[2 * i];
    edges2[2 * i + 1] = edges[2 * i + 1];
    weight2[i] = weight[i];
    indices[i] = i;
  }

  /* While there is a bad edge in G2, remove it                               */

  bad_edge = DPfindBadEdgeK (nnodes, *nedges2, edges2, weight2, k_param);

  while (bad_edge != -1)
  {

    fprintf (stdout, "SEE: Remove edge (%6d %6d) [w = %lf]\n",
             edges2[2 * bad_edge], edges2[2 * bad_edge + 1], weight2[bad_edge]);
    fflush (stdout);

    /*      
     * bad_edge_2 = DPfindBadBinEdgeK(nnodes, *nedges2, edges2, weight2, k_param);
     * 
     * fprintf(stdout, "SEE: BIN   says: remove edge (%6d %6d) [w = %lf] [i = %d]\n", 
     * edges2[2*bad_edge_2], edges2[2*bad_edge_2+1], weight2[bad_edge_2], bad_edge_2);
     * fflush(stdout);
     */

    elim_indices[*nelim] = indices[bad_edge];
    (*nelim) += 1;

    iswap = indices[bad_edge];
    indices[bad_edge] = indices[*nedges2 - 1];
    indices[*nedges2 - 1] = iswap;

    iswap = edges2[2 * bad_edge];
    edges2[2 * bad_edge] = edges2[2 * (*nedges2 - 1)];
    edges2[2 * (*nedges2 - 1)] = iswap;

    iswap = edges2[2 * bad_edge + 1];
    edges2[2 * bad_edge + 1] = edges2[2 * (*nedges2 - 1) + 1];
    edges2[2 * (*nedges2 - 1) + 1] = iswap;

    dswap = weight2[bad_edge];
    weight2[bad_edge] = weight2[*nedges2 - 1];
    weight2[*nedges2 - 1] = dswap;

    (*nedges2) -= 1;

    bad_edge = DPfindBadEdgeK (nnodes, *nedges2, edges2, weight2, k_param);

  }

  free (indices);

#if DISPLAY_MODE
  fprintf (stdout, "\rBIN: Eliminated %d edges.\n", *nelim);
  fprintf (stdout, "BIN: Completed planar sub-graph search.\n");
  fflush (stdout);
#endif

  return 0;

}

int DPfindBadBinEdge (int nnodes,
                      int nedges,
                      int *edges)
{

  /* local variables */
  EGmemPool_t *mem;
  int ltop,
    lbottom,
    lmiddle;
  int is_G_planar_b,
    res = INT_MAX;

  /* basic set-up */
  mem = EGnewMemPool (8192, EGmemPoolNewSize, 4096);
  is_G_planar_b = isPlanarBoyer (nnodes, nedges, edges);
  if (is_G_planar_b)
    goto CLEANUP;

  /* prepare to the binary search */
  ltop = nedges - 1;
  lbottom = 0;
  lmiddle = lbottom + (ltop - lbottom) / 2;

  while (ltop != lbottom)
  {

    /* check if adding it makes the graph non-planar */
    is_G_planar_b = isPlanarBoyer (nnodes, lmiddle + 1, edges);

    /* if so, leave them and move lbottom up. */
    if (is_G_planar_b)
      lbottom = lmiddle + 1;

    /* if not, remove them and move ltop down */
    else
      ltop = lmiddle;

    lmiddle = lbottom + (ltop - lbottom) / 2;

  }                             /* end while */

  is_G_planar_b = isPlanarBoyer (nnodes, lmiddle + 1, edges);

  /* We found the edge that needs to be eliminated */
  if (lbottom < nedges)
    res = lbottom;

  /* ending */
CLEANUP:
  EGfreeMemPool (mem);
  return res;
}

int DPfindBadBinEdgeK (int nnodes,
                       int nedges,
                       int *edges,
                       double *weight,
                       int k)
{

  int rval,
    bad_edge = -1;
  int who,
    nmedges = 0,
   *medges;
  double rvalue;

  if (k == 1)
    return (DPfindBadEdge (nnodes, nedges, edges, weight));
  else if (k < 1)
    return -1;

  medges = (int *) malloc (sizeof (int) * nedges);

  rval = isPlanarBoyer (nnodes, nedges, edges);

  if (rval)
    goto CLEANUP;

  if (k > nedges)
    k = nedges;

  //rval = DPgetBinMinor(nnodes, nedges, edges, &nmedges, medges, nedges);
  rval = DPgetBinMinor (nnodes, nedges, edges, &nmedges, medges, k);
  CHECKRVAL (rval);

  //fprintf(stdout, "DPfindBadBinEdgeK: nmedges = %d\n", nmedges);

  rvalue = random ();
  rvalue /= EGRAND_MAX;

  who = (int) (k * rvalue);

  if (who == k)
    who = k - 1;

  bad_edge = medges[who];

CLEANUP:

  free (medges);

  return bad_edge;

}

int DPgetBinMinor (int nnodes,
                   int nedges,
                   int *edges,
                   int *nmedges,
                   int *medges,
                   int k)
{

  /* local variables */

  EGmemPool_t *mem;
  mem = EGnewMemPool (512, EGmemPoolNewSize, 4096);

  int badedge = 0,
   *ledges,
   *lind,
    lnedges = nedges,
    itmp;
  register int i;

  /* basic set-up */
  *nmedges = 0;
  ledges = (int *) malloc (sizeof (int) * nedges * 2);
  lind = (int *) malloc (sizeof (int) * nedges);
  memcpy (ledges, edges, sizeof (int) * nedges * 2);
  for (i = nedges; i--;)
    lind[i] = i;

  while ((*nmedges < lnedges) && (*nmedges < k))
  {

    /* get next */
    badedge = DPfindBadBinEdge (nnodes, lnedges, ledges);
    if (badedge >= lnedges || badedge < *nmedges)
      break;

    /* swap end 0 */
    itmp = ledges[badedge * 2];
    ledges[badedge * 2] = ledges[(*nmedges) * 2];
    ledges[(*nmedges) * 2] = itmp;

    /* swap end 1 */
    itmp = ledges[badedge * 2 + 1];
    ledges[badedge * 2 + 1] = ledges[(*nmedges) * 2 + 1];
    ledges[(*nmedges) * 2 + 1] = itmp;

    /* save index */
    medges[(*nmedges)] = lind[badedge];

    /* swap index */
    itmp = lind[*nmedges];
    lind[*nmedges] = lind[badedge];
    lind[badedge] = itmp;

    /* increment */
    (*nmedges) += 1;
    lnedges = badedge + 1;

  }                             /* end while */

  /* ending */
  EGfreeMemPool (mem);
  free (ledges);
  free (lind);
  return 0;
}

---