eg_dpseparator.c

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>

#include "eg_list.h"
#include "eg_mempool.h"
#include "graph_boyer.h"
#include "graphdual.h"

#include "eg_util.h"
#include "eg_dgraph.h"
#include "eg_ddomino.h"
#include "eg_pdp.h"
#include "eg_util.h"
#include "eg_ddpconstraint.h"
#include "bc_spanning.h"
#include "bc_util.h"
#include "cookInterface.h"
#include "eg_1pchecker.h"
#include "eg_emptyhandles.h"
double DDP_HEURISTIC_MAXTIME = 10.0;

/* Function Definitions */
int DPseparator (int nnodes,
                 int norig_edges,
                 int *const orig_edges,
                 double *const orig_weight,
                 int *const nIneq,
                 int **const ndominoes,
                 int ***const naset,
                 int ***const nbset,
                 int **const nhandle,
                 int ****const aset,
                 int ****const bset,
                 int ***const handle,
                 const char *const boss_name,
                 double percentage,
                 double ddp_heuristic_maxtime)
{
  /* auxiliary variables */
  int i,
    rval;
  int is_connected;
  /* flags */
  int eliminated_edges_b = 0,
    is_G_planar_b = 0;
  int l_count = 0;
  unsigned int n_empty_handle = 0;
  /* where we will store the edges that are bigger than ZERO_X_EPSILON */
  int *edges = 0,
    nedges = 0,
    nelim_edges = 0;
  double *weight = 0;
  double dpViolation;
  double *slack = 0;
  /* for the primal-dual-primal conversions we need an 'applegate' graph */
  graphP primalG = 0;
  /* the bi-directed dual over which we will search for ddominoes */
  EGdGraph_t *biDualG = 0;
  /* the list of ddominoes */
  EGlist_t *dlist = 0;
  EGlist_t **dembed = 0;
  /* the list of dual dp constraints */
  EGlist_t *ddpc_list = 0;
  /* memory pool */
  EGmemPool_t *mem = 0;
  /* used for handling the list of ddominoes */
  EGddomino_t *ddom;
  EGlistNode_t *ddom_it,
   *ddp_it;
#if ELIM_EDGES
  int *planar_edges = 0,
    nplanar_edges = 0,
   *elim_edges = 0;
  double *planar_weight = 0,
    error = 0;
#endif
#if DISPLAY_MODE
  int nfeasible = 0;
  int ncorrected = 0;
  int nnegative = 0;
  double dual_val,
    prim_val;
#endif
  /* seting up time */
  DDP_HEURISTIC_MAXTIME = ddp_heuristic_maxtime;
  fprintf (stderr, "DP Version %s\n", DP_VERSION_STRING);
#if DISPLAY_MODE
  if (!boss_name)
    fprintf (stdout, "DP Separator: NO boss.\n");
  else
    fprintf (stdout, "DP Separator: Noss name = %s.\n", boss_name);
  fflush (stdout);
#endif

#if LOG_MODE
  saveBgraph ("graph.b.x", nnodes, norig_edges, orig_edges, orig_weight);
#endif

  /* initialize memory pool                                                   */
  mem = EGnewMemPool (512, EGmemPoolNewSize, EGmemPoolNewSize (1));

  /* initialize linked lists                                                  */
  dlist = EGnewList (mem);
  ddpc_list = EGnewList (mem);

#if DISPLAY_MODE
  fprintf (stdout, "DP Separator: Graph G^* has %d nodes, and %d edges.\n",
           nnodes, norig_edges);
  fprintf (stdout,
           "DP Separator: Number of necessary edges for non-planarity = %d.\n",
           3 * nnodes - 5);
  fflush (stdout);
#endif

  /* remove small edges (ie: those smaller than ZERO_X_EPSILON)               */
  rval = removeSmallEdges (ZERO_X_EPSILON, orig_edges, norig_edges, orig_weight,
                           &nedges, &edges, &weight, mem);
  CHECKRVAL (rval);

#if DISPLAY_MODE
  fprintf (stdout,
           "DP Separator: Eliminated %d edges of value less than %lf.\n",
           (norig_edges - nedges), ZERO_X_EPSILON);
  fflush (stdout);
  fprintf (stdout, "DP Separator: Graph G^* now has %d nodes, and %d edges.\n",
           nnodes, nedges);
  fflush (stdout);
#endif

  /* check if the graph is planar                                             */
  primalG = gp_New ();
  rval = cook2boyerGraph (primalG, nnodes, nedges, edges);
  CHECKRVAL (rval);
  is_G_planar_b = gp_Embed (primalG, EMBEDFLAGS_PLANAR);
  rval = gp_SortVertices (primalG);
  EXITRVAL (rval);

  /* if not planar, run heuristic to planarize                                */
  if (is_G_planar_b != OK)
  {
#if DISPLAY_MODE
    fprintf (stdout, "DP Separator: Graph G^* is not planar.\n");
    fflush (stdout);
#endif
    /* First: Contractions.                                                   */
    /* Second: Eliminate Edges using the Kruskal Heuristic.                   */
#if ELIM_EDGES
#if DISPLAY_MODE
    fprintf (stdout, "DP Separator: Running edge-elimination "
             "heuristic to find planar sub-graph.\n");
    fflush (stdout);
#endif

    /* cuidado! */
    planar_edges = EGmemPoolMalloc (mem, sizeof (int) * 2 * nedges);
    planar_weight = EGmemPoolMalloc (mem, sizeof (double) * nedges);
    elim_edges = EGmemPoolMalloc (mem, sizeof (int) * nedges);
    /* cuidado! */

    rval =
      DPedgeEliminationHeuristic (nnodes, nedges, edges, weight, &nplanar_edges,
                                  planar_edges, planar_weight, &nelim_edges,
                                  elim_edges);

    /*  
     * rval = 
     * DPfastEdgeEliminationHeuristic(nnodes, nedges, edges, weight, &nplanar_edges,
     * planar_edges, planar_weight, &nelim_edges, elim_edges, 10);
     */

    CHECKRVAL (rval);

#if DISPLAY_MODE
    fprintf (stdout, "DP Separator: Completed edge-elimination.");
    fflush (stdout);
#endif

    eliminated_edges_b = 1;
    TEST (nplanar_edges == nedges, "ERROR: Graph is not planar, yet no edges "
          "were eliminated.");

    for (i = 0; i < nedges - nplanar_edges; i++)
      error += weight[elim_edges[i]];

#if DISPLAY_MODE
    fprintf (stdout, "DP Separator: Found a planar sub-graph with %d edges.\n",
             nplanar_edges);
    fflush (stdout);
    fprintf (stdout, "DP Separator: Total weight of eliminated edges is %lf.\n",
             error);
    fflush (stdout);
#endif
#endif
  }

  /* if graph is planar, we use the original edge set.                        */
  else
  {
#if DISPLAY_MODE
    fprintf (stdout, "DP Separator: Graph G^* is planar.\n");
    fflush (stdout);
#endif
    eliminated_edges_b = 0;
    planar_edges = edges;
    nplanar_edges = nedges;
    planar_weight = weight;
  }

#if LOG_MODE
  saveBgraph ("graph.planar.b.x", nnodes, nplanar_edges, planar_edges,
              planar_weight);
  saveGraph ("graph.planar.x", nnodes, nplanar_edges, planar_edges,
             planar_weight);
#endif
#if DISPLAY_MODE
  fprintf (stdout, "DP Separator: Computing dual of G^*.\n");
  fflush (stdout);
#endif
  is_connected = EGisCookGraphConnected (nnodes, nplanar_edges, planar_edges,
                                         mem);

#if DISPLAY_MODE
  if (is_connected)
    fprintf (stdout, "DP Separator: G^* is connected.\n");
  else
    fprintf (stdout, "DP Separator: G^* is not connected.\n");
#endif

  if (!is_connected)
    goto END_DP;

  /* compute dual of graph                                                    */
  gp_Free (&primalG);
  primalG = gp_New ();
  rval = cook2boyerGraph (primalG, nnodes, nplanar_edges, planar_edges);
  CHECKRVAL (rval);
  biDualG = getDualBoyer (primalG, nnodes, nplanar_edges, planar_weight,
                          &dembed, mem,0);
  if (!biDualG)
    goto END_DP;

#if DISPLAY_MODE
  fprintf (stdout, "DP Separator: Bi-directed dual graph "
           "has %d nodes and %d edges.\n",
           biDualG->nodes->size, biDualG->nedges);
  fprintf (stdout, "DP Separator: Looking for dual dominoes.\n");
  fflush (stdout);
#endif

  /* with dual of planar graph, obtain dual dominoes                          */
  if (boss_name == 0)
    rval = EGddominoComputeAll (mem, biDualG, dlist, 1, percentage);
  else
    rval = EGddominoComputeAllRemote (mem, biDualG, dlist, 1, boss_name,
                                      DDP_HEURISTIC_MAXTIME);
  CHECKRVAL (rval);

#if DISPLAY_MODE
  fprintf (stdout, "DP Separator: Found %d dual dominoes.\n", dlist->size);
  fflush (stdout);
#endif

  /* now we untengle all the obtained dominoes */
  rval = EGuntangleAllDomino (dlist, biDualG, dembed, mem);
  CHECKRVAL (rval);

  if (!dlist->size)
    goto END_DP;

  /* if using an edge-elimination heuristic, we might want to fix (some? all?) domino values and store this corrected value in primalValue. If there is no edge-elimination we also need to set primalValue, however, in this case there is no need for corrections. Finally, we need to fix numerical errors. Specifically, should we encounter dominoes with slightly negative values, we will round these up to zero (slightly negative == DP_RESET_VALUE. */
  rval =
    EGfixDualDominos (dlist, nnodes, nedges, weight, edges,
                      elim_edges, nelim_edges, primalG, mem);
  CHECKRVAL (rval);


#if DISPLAY_MODE
  if (eliminated_edges_b)
  {

    ncorrected = 0;
    nnegative = 0;
    nfeasible = 0;

    for (ddom_it = dlist->begin; ddom_it; ddom_it = ddom_it->next)
    {

      ddom = (EGddomino_t *) (ddom_it->this);
      dual_val = EGdijkstraCostToLf (ddom->value);
      prim_val = EGdijkstraCostToLf (ddom->primalValue);

      if (fabs (dual_val - prim_val) > 0.0001)
      {
        ncorrected++;
      }
      if (dual_val < 0)
      {
        nnegative++;
      }
      if (prim_val < 1.0000)
        nfeasible++;
    }
    fprintf (stdout,
             "DP Separator: Corrected %d values. %d / %d dominoes remain feasible.\n",
             ncorrected, nfeasible, dlist->size);
    fflush (stdout);
  }
#endif

  /* Solve the cycle problems to compute the dual form of the dp constraints  */
  rval = EGddpcComputeAll (mem, biDualG, dlist, ddpc_list, nedges, 1);
  CHECKRVAL (rval);

#if DISPLAY_MODE

  double thresh,
    thresh_n;
  unsigned int range[DISPLAY_RESOLUTION],
    max_violated = 0;

  for (i = 0; i < DISPLAY_RESOLUTION; i++)
    range[i] = 0;

#endif

  /* reset all values to NULL */
  *nIneq = 0;
  *ndominoes = 0;
  *naset = 0;
  *nbset = 0;
  *aset = 0;
  *bset = 0;
  *nhandle = 0;
  *handle = 0;

  /* ask for space for the constraints */
  if (ddpc_list->size)
  {
    *ndominoes = EGsMalloc (int,
                            ddpc_list->size);
    *naset = EGsMalloc (int *,
                        ddpc_list->size);
    *nbset = EGsMalloc (int *,
                        ddpc_list->size);
    *nhandle = EGsMalloc (int,
                          ddpc_list->size);
    *aset = EGsMalloc (int **,
                       ddpc_list->size);
    *bset = EGsMalloc (int **,
                       ddpc_list->size);
    *handle = EGsMalloc (int *,
                         ddpc_list->size);
    (*aset)[0] = 0;
  }

  if (ddpc_list->size)
    slack = EGsMalloc (double,
                       ddpc_list->size);

#if DISPLAY_MODE
  fprintf (stdout,
           "DP Separator: Primalizing constraints and computing primal slack.\n");
#endif

  for (ddp_it = ddpc_list->begin; ddp_it; ddp_it = ddp_it->next)
  {
    rval = EGgetPrimalDP ((EGddpConstraint_t *) (ddp_it->this),
                          nnodes,
                          nedges,
                          edges,
                          elim_edges,
                          nelim_edges,
                          weight,
                          (*ndominoes) + (*nIneq),
                          (*naset) + (*nIneq),
                          (*nbset) + (*nIneq),
                          (*nhandle) + (*nIneq),
                          (*aset) + (*nIneq),
                          (*bset) + (*nIneq),
                          (*handle) + (*nIneq), &dpViolation, primalG, mem);
    CHECKRVAL (rval);
    if (!((*nhandle)[*nIneq]))
      n_empty_handle++;
    //fprintf(stderr,"%p %d::%d::%u no\n",(void*)((*aset)[0]), *nIneq,l_count,ddpc_list->size);
    l_count++;
    if (dpViolation < 0)
    {
      slack[*nIneq] = dpViolation;
      (*nIneq)++;
#if DISPLAY_MODE
      for (i = 0; i < DISPLAY_RESOLUTION; i++)
      {
        thresh = IDR * i + IDR;
        thresh_n = IDR * i;
        if (thresh_n <= -1 * dpViolation && -1 * dpViolation <= thresh)
        {
          range[i] += 1;
          break;
        }
      }

      if (-1 * dpViolation > 1.0 - DP_MAXVIOLATED_EPSILON)
        max_violated += 1;

#endif

    }                           /* end case we have found a violated constraint */
    else
    {
      for (i = (*ndominoes)[*nIneq]; i--;)
      {
        EGfree ((*aset)[*nIneq][i]);
        EGfree ((*bset)[*nIneq][i]);
      }
      EGfree ((*naset)[*nIneq]);
      EGfree ((*nbset)[*nIneq]);
      EGfree ((*aset)[*nIneq]);
      EGfree ((*bset)[*nIneq]);
      if ((*nhandle)[*nIneq])
        EGfree ((*handle)[*nIneq]);
    }                           /* end freing the non violated primal constraint */
  }                             /* end for ddp_it */

#if DISPLAY_MODE
  fprintf (stdout, "DP Separator: Found %d primal DP-Constraints.\n", *nIneq);
  fprintf (stdout, "DP Separator: Found %u empty handles.\n", n_empty_handle);
  fflush (stdout);

  if (nelim_edges || 1)
  {

    fprintf (stdout, "DP Separator: Final |slack values|\n");
    fflush (stdout);

    for (i = 0; i < DISPLAY_RESOLUTION; i++)
    {
      thresh = IDR * i + IDR;
      thresh_n = IDR * i;
      if (range[i])
        fprintf (stdout, "[%8lf, %8lf]   %d\n", thresh_n, thresh, range[i]);
      fflush (stdout);
    }

    fprintf (stdout, "Max-Violations        %d\n", max_violated);
    fflush (stdout);

  }
#endif

  if (*nIneq)
  {
    double *violation;
    int **e_coeff;
    violation = EGsMalloc (double,
                           *nIneq);
    e_coeff = EGsMalloc (int *,
                         *nIneq);
    for (i = *nIneq; i--;)
      e_coeff[i] = EGsMalloc (int,
                              norig_edges);
    rval = EG1pChecker (nnodes, norig_edges, orig_edges, orig_weight, *nIneq,
                        *ndominoes, *naset, *nbset, *nhandle,
                        *aset, *bset, *handle, violation, e_coeff);
    CHECKRVAL (rval);
    for (i = *nIneq; i--;)
    {
      WARNING (fabs (violation[i] - slack[i]) > 1e-2,
               "Violation don't agree for ineq "
               "%i\n(computed)%7.5lf (alg)%7.5lf (diff)%7.5lf\n", i,
               violation[i], slack[i], fabs (violation[i] - slack[i]));
      EGfree (e_coeff[i]);
    }
    EGfree (violation);
    EGfree (e_coeff);
  }
  if (slack)
    EGfree (slack);

#if DP_REMOVE_EMPTY_HANDLES
  if (n_empty_handle)
  {

#if DISPLAY_MODE
    fprintf (stdout, "DP Separator: Removing empty-handle constraints.\n");
    fflush (stdout);
#endif

    rval =
      DPremoveEmptyHandles (nIneq, ndominoes, naset, nbset, nhandle, aset, bset,
                            handle);
    CHECKRVAL (rval);

#if DISPLAY_MODE
    fprintf (stdout,
             "DP Separator: Final number of non-empty-handle constraints: %d.\n",
             *nIneq);
    fflush (stdout);
#endif

  }
#endif

  /* EGfree up used memory */

END_DP:

#if DISPLAY_MODE
  fprintf (stdout, "DP Separator: Bye.\n");
  fflush (stdout);
#endif

  /* EGfree dlist, ddplist, and pdplist */
  if (dlist)
  {
    if (dlist->size)
      EGlistClearMP (dlist, EGfreeDdomino, mem);
    EGfreeList (dlist);
  }

  if (ddpc_list)
  {
    if (ddpc_list->size)
      EGlistClearMP (ddpc_list, EGfreeDDPconstraintMP, mem);
    EGfreeList (ddpc_list);
  }

  /* EGfree 'graphdual' objects */
  if (weight)
    EGmemPoolFree (weight, sizeof (double) * nedges, mem);
  if (edges)
    EGmemPoolFree (edges, sizeof (int) * 2 * nedges, mem);

#if ELIM_EDGES
  /* EGfree planar-graph heuristic data structures */
  if (eliminated_edges_b)
  {
    EGmemPoolSFree (planar_edges, int,
                    2 * nedges,
                    mem);
    EGmemPoolSFree (elim_edges, int,
                    nedges,
                    mem);
    EGmemPoolSFree (planar_weight, double,
                    nedges,
                    mem);
  }
#endif

  /* EGfree biDualG: weight, edges, and graph structure. */
  if (biDualG)
  {
    EGdGraphClearMP (biDualG, EGfreeMengerEdgeDataMP, EGfreeMengerNodeDataMP, 0,
                     mem, mem, 0);
    EGfreeDGraph (biDualG);
  }

  if (dembed)
  {
    for (i = nplanar_edges - nnodes + 2; i--;)
      EGfreeList (dembed[i]);
    EGmemPoolSFree (dembed, EGlist_t *, nplanar_edges - nnodes + 2, mem);
  }
  /* EGfree memory pool */
  EGfreeMemPool (mem);
  gp_Free (&primalG);

  /* return */
  return 0;

}

---