eg_ddpconstraint.c

Go to the documentation of this file.

#include "eg_ddpconstraint.h"
/* this macro return the 'edge' ponter asociated to the 'b'-th arc in path 'a'
 * of the domino 'dom'. It takes care of different types of dominos in an easy
 * form */
#define EGddGetEdge(DOM,a,b) ((graphNodeP)(((DOM)->DDtype == DOM_CC_ONE) ? \
    (DOM)->path[a][b]:((EGmengerEdgeData_t*)(((EGdGraphEdge_t*)\
    ((DOM)->path[a][b]))->data))->data))

static inline int EGgetPrimalDPdist (EGddpConstraint_t const *const ddp,  /* dual DP constraint */

                                     double *const dist)  /* where to store the distance to the
                                                           * fractional X */
{

  /* local variables */
  double alphax = -1.0;
  register unsigned int k;


  /* now we update the parity of the edges in F */
  for (k = ddp->nF; k--;)
    alphax +=
      EGdijkstraCostToLf (((EGmengerEdgeData_t *) (ddp->F[k]->data))->cost);

  /* here we compute the primal dominoes, and loop through all of them */
  for (k = ddp->ndom; k--;)
    alphax += EGdijkstraCostToLf (ddp->ddom[k]->value);

  alphax /= sqrt (((double) (ddp->nF + ddp->ndom)));
  *dist = alphax;
  return 0;
}

static int EGdominoCompare (const void *p1,
                            const void *p2)
{

  if ((*((ccddomp_t) (p1)))->id < (*((ccddomp_t) (p2)))->id)
    return -1;
  else if ((*((ccddomp_t) (p1)))->id > (*((ccddomp_t) (p2)))->id)
    return 1;
  return 0;
}

static inline int EGgetPrimalDPangle (EGddpConstraint_t const *const ddp1,  /* dual DP constraint */

                                      EGddpConstraint_t const *const ddp2,  /* dual DP constraint */

                                      double *const angle)  /* where to store the distance to the
                                                             * fractional X */
{

  /* we define the norm of a DP inequality as nF + sum ddom_slack^2 */
  /* local variables */
  register unsigned int k1 = 0,
    k2 = 0;
  double norm3 = 0;

  /* now we compute  */
  while (k1 < ddp1->nF && k2 < ddp2->nF)
  {
    switch (EGedgeCompare (ddp1->F + k1, ddp2->F + k2))
    {
    case 0:
      norm3 += 1;
      k1++;
      k2++;
      break;
    case 1:
      k1++;
      break;
    case -1:
      k2++;
      break;
    }
  }
  k1 = k2 = 0;
  while (k1 < ddp1->ndom && k2 < ddp2->ndom)
  {
    switch (EGdominoCompare (ddp1->ddom + k1, ddp2->ddom + k2))
    {
    case 0:
      norm3 += 1;
      k1++;
      k2++;
      break;
    case 1:
      k1++;
      break;
    case -1:
      k2++;
      break;
    }
  }

  *angle = norm3;
  *angle /= sqrt ((double) ((ddp1->nF + ddp1->ndom) * (ddp2->nF + ddp2->ndom)));

  /* ending */
  return 0;
}

static long long unsigned int global_cnt = 0;

static inline double EGddpHeuristicVal (EGdijkstraCost_t dcost)
{

  double val = EGdijkstraCostToLf (dcost);

  if (val < 0.01)
    val = 0.01;

  /* original values:
   * if (val < 0.01)
   * val = 0.01;
   */

  val = (1 / val);

  return val;

}


/* this data structure is used to store relevant information 
   for the edges in a cycle graph */

static inline EGcycleData_t *EGnewCycleData (EGmemPool_t * mem,
                                             EGdijkstraCost_t new_weight,
                                             EGddomino_t * new_ddom,
                                             EGdGraphEdge_t * new_e)
{

  double val;
  EGcycleData_t *cycle_data;

  if (EGdijkstraCostIsLess (new_weight, EGdijkstraToCost (-0.02)))
  {
    fprintf (stdout,
             "WARNING! cycle-data edge is less than -0.02 (%lf) set to 0.0\n",
             EGdijkstraCostToLf (new_weight));
    new_weight = EGdijkstraToCost (0.0);
  }
  cycle_data = EGmemPoolSMalloc (mem, EGcycleData_t, 1);
  cycle_data->weight = new_weight;
  cycle_data->e = new_e;
  cycle_data->ddom = new_ddom;

  val = round (EGddpHeuristicVal (new_weight));

  cycle_data->pweight = ((unsigned int) val);

  return cycle_data;

}

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

  EGmemPoolSFree (v, EGcycleData_t, 1, mem);

  return;

}

EGdGraph_t *EGnewCycleGraph (EGmemPool_t * mem,
                             EGlist_t * dlist,
                             EGdGraph_t * bdG,
                             EGdijkstraCost_t max_val)
{

  unsigned int p,
    ph,
    pt;
  EGdGraphNode_t *n;
  EGdGraphEdge_t *e;
  EGlistNode_t *n_it,
   *e_it,
   *dd_it;
  EGdGraph_t *cycleG;
  EGcycleData_t *edata;
  EGddomino_t *ddom;
  EGdijkstraNodeData_t *ndata;
  EGmengerEdgeData_t *mdata;

  EGdGraphNode_t **node_array;

  cycleG = EGnewDGraph (mem);

  /* add the nodes: there are two nodes in cycleG per every node in bdG       */
  node_array = EGmemPoolSMalloc (mem, EGdGraphNode_t *, 2 * (bdG->nodes->size));

  for (n_it = bdG->nodes->begin; n_it; n_it = n_it->next)
  {
    n = (EGdGraphNode_t *) (n_it->this);

    /* we make a 'left' copy                                                  */
    p = 2 * (n->id);
    ndata = EGnewDijkstraNodeData (mem);
    ndata->father = 0;
    node_array[p] = EGdGraphAddNode (cycleG, ndata);
    /* and a 'right' copy                                                     */
    p = 2 * (n->id) + 1;
    ndata = EGnewDijkstraNodeData (mem);
    ndata->father = 0;
    node_array[p] = EGdGraphAddNode (cycleG, ndata);
  }

  /* add the F edges: there are two F edges per every edge in bdG             */
  for (n_it = bdG->nodes->begin; n_it; n_it = n_it->next)
  {
    n = (EGdGraphNode_t *) (n_it->this);
    for (e_it = n->out_edges->begin; e_it; e_it = e_it->next)
    {
      e = (EGdGraphEdge_t *) (e_it->this);
      mdata = (EGmengerEdgeData_t *) (e->data);
      if (EGdijkstraCostIsLess (max_val, mdata->cost))
        continue;

      /* we make a 'left copy'                                                */

      ph = 2 * (e->head->id);
      pt = 2 * (e->tail->id);
      edata = EGnewCycleData (mem, mdata->cost, 0, e);
      EGdGraphAddEdge (cycleG, edata, node_array[ph], node_array[pt]);

      /* and a 'right' copy                                                   */

      ph = 2 * (e->head->id) + 1;
      pt = 2 * (e->tail->id) + 1;
      edata = EGnewCycleData (mem, mdata->cost, 0, e);
      EGdGraphAddEdge (cycleG, edata, node_array[ph], node_array[pt]);

    }
  }

  /* add the odd edges                                                        */
  for (dd_it = dlist->begin; dd_it; dd_it = dd_it->next)
  {
    ddom = (EGddomino_t *) (dd_it->this);
    if (EGdijkstraCostIsLess (max_val, ddom->primalValue))
      continue;

    /* we make a 'left s' to 'right t' copy                                   */

    ph = 2 * (ddom->s->id);
    pt = 2 * (ddom->t->id) + 1;
    edata = EGnewCycleData (mem, ddom->primalValue, ddom, 0);
    EGdGraphAddEdge (cycleG, edata, node_array[ph], node_array[pt]);

    /* we make a 'right s' to 'left t' copy                                   */

    ph = 2 * (ddom->s->id) + 1;
    pt = 2 * (ddom->t->id);
    edata = EGnewCycleData (mem, ddom->primalValue, ddom, 0);
    EGdGraphAddEdge (cycleG, edata, node_array[ph], node_array[pt]);

    /* we make a 'left t' to 'right s' copy                                   */

    ph = 2 * (ddom->t->id);
    pt = 2 * (ddom->s->id) + 1;
    edata = EGnewCycleData (mem, ddom->primalValue, ddom, 0);
    EGdGraphAddEdge (cycleG, edata, node_array[ph], node_array[pt]);

    /* we make a 'right t' to 'left s' copy                                   */

    ph = 2 * (ddom->t->id) + 1;
    pt = 2 * (ddom->s->id);
    edata = EGnewCycleData (mem, ddom->primalValue, ddom, 0);
    EGdGraphAddEdge (cycleG, edata, node_array[ph], node_array[pt]);
  }

  EGmemPoolSFree (node_array, EGdGraphNode_t *, 2 * (bdG->nodes->size), mem);

  return cycleG;

}

void EGfreeDDPconstraintMP (void *v,
                            EGmemPool_t * mem)
{
  EGddpConstraint_t *ddpc;

  ddpc = (EGddpConstraint_t *) (v);

  if (ddpc->ndom)
    EGmemPoolSFree (ddpc->ddom, EGddomino_t *, ddpc->ndom, mem);
  if (ddpc->nF)
    EGmemPoolSFree (ddpc->F, EGdGraphEdge_t *, ddpc->nF, mem);

  EGmemPoolSFree (v, EGddpConstraint_t, 1, mem);

  return;
}


int EGboyerEdgeNumber (const EGdGraphEdge_t * e)
{
  EGmengerEdgeData_t *d;
  graphNodeP a;

  d = (EGmengerEdgeData_t *) (e->data);
  a = (graphNodeP) (d->data);

  return (a->id);
}

int EGedgeCompare (const void *p1,
                   const void *p2)
{
  const EGdGraphEdge_t *e1,
   *e2;

  e1 = (*((ccep_t) (p1)));
  e2 = (*((ccep_t) (p2)));

  if (EGboyerEdgeNumber (e1) < EGboyerEdgeNumber (e2))
    return -1;
  else if (EGboyerEdgeNumber (e1) > EGboyerEdgeNumber (e2))
    return 1;
  return 0;
}

/* this macro return the dijksra cost of a dual edge */
#define EGddGetEdgeCost(dedge) (\
    ((EGmengerEdgeData_t*)(((EGdGraphEdge_t*)(dedge))->data))->cost)
static unsigned int __ndisc = 0;
static unsigned int __accept = 0;
static inline int EGddpMuCondition (EGddpConstraint_t * ddpc)
{
  EGdijkstraCost_t muH = EGdijkstraToCost (ddpc->ndom + 1),
    p_cost[2];
  register unsigned int i;
  unsigned pth,
    dom;
  /* now we substract the central path of all dominoes and F */
  for (i = ddpc->nF; i--;)
    muH = EGdijkstraCostSub (muH, EGddGetEdgeCost (ddpc->F[i]));
  for (dom = ddpc->ndom; dom--;)
  {
    p_cost[0] = EGdijkstraToCost (0);
    p_cost[1] = EGdijkstraToCost (100);
    /* compute cheapest path in p_cost[1] */
    for (pth = 3; pth--;)
    {
      for (i = ddpc->ddom[dom]->npath[pth]; i--;)
        p_cost[0] = EGdijkstraCostAdd (p_cost[0],
                                       EGddGetEdgeCost (ddpc->ddom[dom]->
                                                        path[pth][i]));
      if (EGdijkstraCostIsLess (p_cost[0], p_cost[1]))
        p_cost[1] = p_cost[0];
    }                           /* end computing cheapest path */
    muH = EGdijkstraCostSub (muH, p_cost[1]);
  }                             /* end loop through all dominoes */
  /* now we check that muH > 0 */
  if (EGdijkstraCostIsLess (EGdijkstraToCost (0), muH))
  {
    __accept++;
    return 1;
  }
  __ndisc++;
  //fprintf(stderr,"Discarded Handle %.3lf%%\r",  
  //        ((double)(__ndisc))/(__ndisc+__accept) );
  return 0;
}

static inline int EGddpIsMinimal (EGddpConstraint_t * ddpc,
                                  unsigned int *marray,
                                  unsigned int marray_size)
{

  unsigned int i;

  for (i = 0; i < marray_size; i++)
    marray[i] = 0;

  for (i = 0; i < ddpc->nF; i++)
  {
    marray[ddpc->F[i]->head->id] += 1;
    marray[ddpc->F[i]->tail->id] += 1;
  }

  for (i = 0; i < ddpc->ndom; i++)
  {
    marray[ddpc->ddom[i]->s->id] += 1;
    marray[ddpc->ddom[i]->t->id] += 1;
  }

  for (i = 0; i < marray_size; i++)
  {
    if (marray[i] > 2)
      return 0;
  }

  return 1;

}

static inline int EGddpExtractSolution (EGdGraphEdge_t ** path,
                                        unsigned int *npath,
                                        EGdGraphNode_t * t,
                                        EGdijkstraCost_t * LHS_val)
{

  EGdijkstraCost_t dcVal;

  EGdGraphEdge_t *e;
  EGdGraphNode_t *n;
  EGdijkstraNodeData_t *dndata;
  EGcycleData_t *cdata;

  unsigned int pos;

  n = t;
  pos = 0;
  dcVal = EGdijkstraToCost (0.0);

  while (((n->id) / 2 != (t->id) / 2) || !pos)
  {
    dndata = (EGdijkstraNodeData_t *) (n->data);
    e = dndata->father;

    cdata = (EGcycleData_t *) (e->data);
    dcVal = EGdijkstraCostAdd (dcVal, cdata->weight);

    path[pos] = e;
    n = e->tail;

    pos += 1;
  }

  *npath = pos;
  *LHS_val = dcVal;

  return 0;

}

static inline int EGddpAreSame (EGddpConstraint_t * c1,
                                EGddpConstraint_t * c2)
{
  unsigned int i;

  if ((c1->nF != c2->nF) || (c1->ndom != c2->ndom))
    return 0;

  for (i = 0; i < c1->nF; i++)
  {
    if (EGboyerEdgeNumber (c1->F[i]) != EGboyerEdgeNumber (c2->F[i]))
      return 0;
  }
  for (i = 0; i < c1->ndom; i++)
  {
    if ((c1->ddom[i]->s->id) != (c2->ddom[i]->s->id))
      return 0;
    if ((c1->ddom[i]->t->id) != (c2->ddom[i]->t->id))
      return 0;
    if ((c1->ddom[i]->value != c2->ddom[i]->value))
      return 0;
  }

  return 1;

}

static inline int EGddpSortData (EGddpConstraint_t * c)
{
  if (c->ndom > 1)
    qsort (c->ddom, c->ndom, sizeof (EGddomino_t *), EGdominoCompare);
  if (c->nF > 1)
    qsort (c->F, c->nF, sizeof (EGdGraphEdge_t *), EGedgeCompare);

  return 0;
}

static inline EGddpConstraint_t *EGnewDDPconstraint (EGmemPool_t * mem,
                                                     EGdGraphEdge_t ** dij_sol,
                                                     unsigned int ndij_sol,
                                                     int k,
                                                     EGdijkstraCost_t val)
{

  EGcycleData_t *cdata;
  EGddpConstraint_t *ddpc;
  unsigned int i;

  ddpc = EGmemPoolSMalloc (mem, EGddpConstraint_t, 1);
  ddpc->ndom = 0;
  ddpc->nF = 0;

  for (i = 0; i < ndij_sol; i++)
  {
    cdata = (EGcycleData_t *) (dij_sol[i]->data);
    if (cdata->ddom)
      ddpc->ndom += 1;
    else
      ddpc->nF += 1;
  }

  if (ddpc->nF)
    ddpc->F = EGmemPoolSMalloc (mem, EGdGraphEdge_t *, ddpc->nF);
  else
    ddpc->F = 0;
  if (ddpc->ndom)
    ddpc->ddom = EGmemPoolSMalloc (mem, EGddomino_t *, ddpc->ndom);
  else
    ddpc->ddom = 0;

  ddpc->ndom = 0;
  ddpc->nF = 0;

  for (i = 0; i < ndij_sol; i++)
  {
    cdata = (EGcycleData_t *) (dij_sol[i]->data);
    if (cdata->ddom)
    {
      ddpc->ddom[ddpc->ndom] = cdata->ddom;
      ddpc->ndom += 1;
    }
    else
    {
      ddpc->F[ddpc->nF] = cdata->e;
      ddpc->nF += 1;
    }
  }

  ddpc->LHS_dp = val;
  ddpc->slack_dp = EGdijkstraCostToLf (val) - 1.0;
  return ddpc;

}


static int EGddpAddCutIfNew (EGddpConstraint_t ** ddpc_array,
                             double *const ddpc_dist,
                             double *const ddpc_angle_norm,
                             double *const *const ddpc_angle,
                             unsigned int *const ddpc_size,
                             unsigned int const ddpc_max_size,
                             double *const ddpc_best_angle,
                             double *const ddpc_worst_angle,
                             unsigned int *const ddpc_worst_angle_id,
                             EGdGraphNode_t * nrigh,
                             EGdGraph_t * cycleG,
                             EGdGraphEdge_t ** dij_sol,
                             unsigned int *marray,
                             char *sarray,
                             size_t * os,
                             int k,
                             EGmemPool_t * mem)
{

  int rval;
  unsigned int ndij_sol;
  unsigned register int i;
  unsigned int hcmax = 0,
    worst_dist = 0,
    pos = *ddpc_size;
  EGddpConstraint_t *ddpc;
  EGdijkstraCost_t lhs_val;
  double val,
    dist,
    maxangle,
    dtmp,
    norm = 1.0;

  EGddpExtractSolution (dij_sol, &ndij_sol, nrigh, &lhs_val);

  os = 0;
  /* ensure that the constraint is sufficiently violated */
  if (((1.0) - EGdijkstraCostToLf (lhs_val)) <= ((DDP_MIN_VIOLATION)))
    return 0;

  /* create the new ddpc constraint */
  ddpc = EGnewDDPconstraint (mem, dij_sol, ndij_sol, 1, lhs_val);

  /* if the constraint don't have an odd number of doinos, we discardit and
   * return (no error) */
  if (!(ddpc->ndom & 1))
  {
    EGfreeDDPconstraintMP (ddpc, mem);
    return 0;
  }

  /* ensure the constraint got at least 3 dominoes */
  if (ddpc->ndom < 3)
  {
    EGfreeDDPconstraintMP (ddpc, mem);
    return 0;
  }

  rval = EGddpIsMinimal (ddpc, marray, cycleG->nodes->size / 2);
  if (!rval)
  {
    sarray[(nrigh->id) / 2] = 'M';
    EGfreeDDPconstraintMP (ddpc, mem);
    return 0;
  }

  /* sort the internal data in the new constraint */
  rval = EGddpSortData (ddpc);
  CHECKRVAL (rval);

  /* check if the inequality is already added to the pool */
  for (i = *ddpc_size; i--;)
  {
    if (EGddpAreSame (ddpc, ddpc_array[i]))
    {
      sarray[(nrigh->id) / 2] = 'R';
      EGfreeDDPconstraintMP (ddpc, mem);
      return 0;
    }
  }

  /* compute the distance from the fractional point to the newly found
   * constraint */
  EGgetPrimalDPdist (ddpc, &dist);
  val = -1.0 * dist;

  if (pos < ddpc_max_size)
  {
    ddpc_array[pos] = ddpc;
    ddpc_dist[pos] = val;
    ddpc_angle_norm[pos] = 1.0;
    /* compute the angles */
    if (pos)
    {
      for (i = pos; i--;)
      {
        EGgetPrimalDPangle (ddpc, ddpc_array[i], &dtmp);
        ddpc_angle[pos][i] = dtmp;
        ddpc_angle[i][pos] = dtmp;
        ddpc_angle_norm[i] += dtmp * dtmp;
        ddpc_angle_norm[pos] += dtmp * dtmp;
      }
    }
    else
      *ddpc_worst_angle_id = 0;
    (*ddpc_size)++;
    pos++;
    MESSAGE (EG_DP_SELECT_DEBUG, "ddpc added because pool not full");
  }
  /* if the pool is full, then we have to choose who to kick out of the pool,
   * we do that taking into account the angle and the distance */
  else
  {
    maxangle = 0.0;
    worst_dist = ddpc_max_size - 1;
    /* compute  the angle from the new cut to all cuts in the pool */
    for (i = ddpc_max_size; i--;)
    {
      EGgetPrimalDPangle (ddpc, ddpc_array[i], &dtmp);
      ddpc_angle[ddpc_max_size][i] = dtmp;
      norm += dtmp * dtmp;
      if (dtmp > maxangle)
      {
        hcmax = i;
        maxangle = dtmp;
      }
      if (ddpc_dist[worst_dist] > ddpc_dist[i])
        worst_dist = i;
    }
    /* now we  have the closest angle constraint in the pool, we first check if
     * the current cut is near an added cut and its distancce to X is better,
     * if so, we replace the old cut by the newly found cut */
#if DDPC_REPLACE_CLOSE
    if ((maxangle > DP_BIG_ANGLE) && (ddpc_dist[hcmax] < val))
    {
      MESSAGE (EG_DP_SELECT_DEBUG, "ddpc added because dominate another"
               "inequality, angle(%lf) new disr(%lf) old dist(%lf)", maxangle,
               val, ddpc_dist[hcmax]);
      EGfreeDDPconstraintMP (ddpc_array[hcmax], mem);
      ddpc_array[hcmax] = ddpc;
      ddpc_dist[hcmax] = val;
      norm -= ddpc_angle[ddpc_max_size][hcmax] *
        ddpc_angle[ddpc_max_size][hcmax];
      for (i = ddpc_max_size; i--;)
      {
        dtmp = ddpc_angle[i][hcmax];
        ddpc_angle_norm[i] -= dtmp * dtmp;
        dtmp = ddpc_angle[ddpc_max_size][i];
        ddpc_angle_norm[i] += dtmp * dtmp;
        ddpc_angle[hcmax][i] = dtmp;
        ddpc_angle[i][hcmax] = dtmp;
      }
      ddpc_angle_norm[hcmax] = norm;
    }
#else
    if (0)
    {
    }
#endif
    /* if the worst norm is good enough, we put in the inequality if it improves
     * the wordt cuality */
#if DDPC_BEST_VIOLATION
    else if (val > ddpc_dist[worst_dist])
    {
      pos = worst_dist;
      MESSAGE (EG_DP_SELECT_DEBUG, "ddpc added because improve worst distance"
               " with norm(%lf) angle(%lf) new dist(%lf) old dist(%lf)",
               norm, maxangle, val, ddpc_dist[pos]);
      EGfreeDDPconstraintMP (ddpc_array[pos], mem);
      ddpc_array[pos] = ddpc;
      ddpc_dist[pos] = val;
      norm -= ddpc_angle[ddpc_max_size][pos] * ddpc_angle[ddpc_max_size][pos];
      for (i = ddpc_max_size; i--;)
      {
        dtmp = ddpc_angle[i][pos];
        ddpc_angle_norm[i] -= dtmp * dtmp;
        dtmp = ddpc_angle[ddpc_max_size][i];
        ddpc_angle_norm[i] += dtmp * dtmp;
        ddpc_angle[i][pos] = dtmp;
        ddpc_angle[pos][i] = dtmp;
      }
      ddpc_angle_norm[pos] = norm;
    }
#endif
    /* if the cut don't dominate another cut, it still may improve the spread
     * of the set of constraints in the pool, if so, we use it */
#if DDPC_IMPROVE_SPREAD
    else if (*ddpc_worst_angle > SELECT_MIN_NORM && *ddpc_worst_angle +
             (ddpc_angle[ddpc_max_size][*ddpc_worst_angle_id]) *
             (ddpc_angle[ddpc_max_size][*ddpc_worst_angle_id]) > norm &&
             val > DP_MAXDECREASE * ddpc_dist[worst_dist])
    {
      pos = *ddpc_worst_angle_id;
      MESSAGE (EG_DP_SELECT_DEBUG, "ddpc added because improve worst norm(%lf)"
               " with norm(%lf) angle(%lf) new dist(%lf) old dist(%lf)",
               *ddpc_worst_angle, norm, maxangle, val, ddpc_dist[pos]);
      EGfreeDDPconstraintMP (ddpc_array[pos], mem);
      ddpc_array[pos] = ddpc;
      ddpc_dist[pos] = val;
      norm -= ddpc_angle[ddpc_max_size][pos] * ddpc_angle[ddpc_max_size][pos];
      for (i = ddpc_max_size; i--;)
      {
        dtmp = ddpc_angle[i][pos];
        ddpc_angle_norm[i] -= dtmp * dtmp;
        dtmp = ddpc_angle[ddpc_max_size][i];
        ddpc_angle_norm[i] += dtmp * dtmp;
        ddpc_angle[i][pos] = dtmp;
        ddpc_angle[pos][i] = dtmp;
      }
      ddpc_angle_norm[pos] = norm;
    }
#endif
    /* now we just discard */
    else
    {
      MESSAGE (EG_DP_SELECT_DEBUG, "ddpc discarded. worst norm(%lf) with norm"
               "(%lf) angle(%lf) new dist(%lf) old dist(%lf)",
               *ddpc_worst_angle, norm, maxangle, val, ddpc_dist[worst_dist]);
      EGfreeDDPconstraintMP (ddpc, mem);
    }
  }
  /* now we recompute ddpc_worst/best_angle/_id */
  *ddpc_worst_angle = 0.0;
  *ddpc_best_angle = ddpc_max_size;
  if (*ddpc_size)
    for (i = *ddpc_size; i--;)
    {
      dtmp = ddpc_angle_norm[i];
      if (dtmp > *ddpc_worst_angle)
      {
        *ddpc_worst_angle = dtmp;
        *ddpc_worst_angle_id = i;
      }
      if (dtmp < *ddpc_best_angle)
        *ddpc_best_angle = dtmp;
    }

  return 0;
}

static unsigned char *__edge_valid = 0;
static unsigned int __edge_valid_size = 0;
static inline EGdGraphEdge_t *EGddpHeuristicSample (EGdGraphNode_t * n,
                                                    unsigned int
                                                    *marked_node_array,
                                                    unsigned int
                                                    n_marked_dominoes,
                                                    unsigned int *odd_num,
                                                    int *ddom_markers,
                                                    int k)
{

  int is_e_dom,
    head_pid,
    num_feasible = 0,
    dom_min;
  double rvalue;
  unsigned int irvalue,
    max_weight = 0,
    current_weight = 0;
  EGdGraphNode_t *head;
  EGdGraphEdge_t *e = 0,
   *next_e = 0,
   *saved_e = 0;
  EGcycleData_t *cd;
  EGlistNode_t *it;
  /* these two are used to store which edges are feasibles, so we don't work
   * twice */
  unsigned int edge_valid_it;

  /* basic set-up */
  if (__edge_valid_size < n->out_edges->size)
  {
    if(__edge_valid) EGfree (__edge_valid);
    __edge_valid_size = n->out_edges->size;
    __edge_valid = EGsMalloc (unsigned char,
                              __edge_valid_size);
  }
  memset (__edge_valid, 0, sizeof (unsigned char) * n->out_edges->size);

  /* iterate through the edges to compute max_weight */
  for (edge_valid_it = 0, it = n->out_edges->begin; it;
       it = it->next, edge_valid_it++)
  {

    /* compute information relevant to this edge */
    e = (EGdGraphEdge_t *) (it->this);
    cd = (EGcycleData_t *) (e->data);
    head = e->head;
    head_pid = ((head->id) / 2);

    /* check if the edge is forbidden */
    if (cd->ddom && ddom_markers && ddom_markers[cd->ddom->id])
      continue;

    /* check if the edge leads to a previously visited node, 
     * and if so, check if the number of dominoes visited 
     * since then is odd and bigger-than-or-equal to DP_HEUR_MIN_DOM  */
    if (marked_node_array[head_pid])
    {
      rvalue = random ();
      rvalue /= EGRAND_MAX;
      is_e_dom = (cd->ddom ? 1 : 0);
      if (rvalue < DP_HEUR_MORE_DOM_BIAS)
        dom_min = 2 + n_marked_dominoes + is_e_dom - odd_num[head_pid];
      else
        dom_min = n_marked_dominoes + is_e_dom - odd_num[head_pid];
      /* check if this edge would give a valid cycle, of so, with some
       * probability return this edge, otherwise, discard this edge */
      if ((dom_min >= DP_HEUR_MIN_DOM) && (dom_min & 1))
      {
        rvalue = random ();
        rvalue /= EGRAND_MAX;
        if (rvalue < DP_HEUR_CLOSE_BIAS)
        {
          next_e = e;
          goto FOUND_EDGE;
        }
      }
      else
        continue;
    }

    /* we have a feasible edge: update data */
    __edge_valid[edge_valid_it] = 1;
    num_feasible += 1;
    max_weight += cd->pweight;
    saved_e = e;
  }

  /* check the trivial cases */
  if (!num_feasible)
    return 0;

  if (num_feasible == 1)
  {
    next_e = saved_e;
    goto FOUND_EDGE;
  }

  /* sample a value between 0 and max_weight */
  rvalue = random ();
  rvalue /= EGRAND_MAX;
  rvalue *= max_weight;
  rvalue = round (rvalue);
  irvalue = (unsigned int) (rvalue);

  /* iterate through the edges to choose the edge */
  for (edge_valid_it = 0, it = n->out_edges->begin; it;
       it = it->next, edge_valid_it++)
  {
    /* discard non-valid edges */
    if (!__edge_valid[edge_valid_it])
      continue;

    /* compute information relevant to this edge */
    e = (EGdGraphEdge_t *) (it->this);
    current_weight += ((EGcycleData_t *) (e->data))->pweight;
    if (current_weight >= irvalue)
    {
      next_e = e;
      break;
    }

  }

FOUND_EDGE:

  /* make sure we found an edge */
  ADVTESTL (1, !next_e, 0, "error sampling edge");

  return next_e;

}


/* EGddpAddHeuristicCuts_3
 marray:  used to mark which nodes have been traversed by the random walk.
 odd_num: used to indicate how many odd edges have been traversed in order 
 to get to each node */
int EGddpAddHeuristicCuts_3 (EGddpConstraint_t ** ddpc_array,
                             double *const ddpc_dist,
                             double *const ddpc_angle_norm,
                             double *const *const ddpc_angle,
                             unsigned int *const ddpc_size,
                             unsigned int const ddpc_max_size,
                             double *const ddpc_best_angle,
                             double *const ddpc_worst_angle,
                             unsigned int *const ddpc_worst_angle_id,
                             int const nedges,
                             EGdGraphNode_t * start_node,
                             EGdGraphEdge_t * start_edge,
                             int k,
                             double ubound,
                             EGdGraph_t * cycleG,
                             EGdGraphEdge_t ** dij_sol,
                             unsigned int *marray,
                             char *sarray,
                             unsigned int *odd_num,
                             unsigned int nddom,
                             int *ddom_markers,
                             size_t * os,
                             double max_node_time,
                             EGmemPool_t * mem)
{

  int rval;
  unsigned int stop = 0;

  EGdijkstraCost_t path_val;
  unsigned int current_pid,
    n_marked_dominoes;

  EGdGraphNode_t *current_node;
  EGdGraphEdge_t *current_edge;

  EGdijkstraNodeData_t *node_data;
  EGcycleData_t *edge_data;

  EGtimer_t t;

  /* ensure that there is a start node */
  if (!start_node)
    start_node = start_edge->tail;

  /* if shortest path already exceeds bound, we abort */
  if (sarray[(start_node)->id / 2] == 'U')
    return 0;

  EGtimerReset (&t);
  EGtimerStart (&t);

  while (!stop)
  {

    EGtimerStop (&t);

    if (t.time > max_node_time)
    {
      stop = 1;
      break;
    }

    EGtimerStart (&t);

    global_cnt += 1;

    /* reset the current path val to zero */
    path_val = EGdijkstraToCost (0.0);

    /* set all node markers to zero as we have not visited any yet */
    memset (marray, 0, sizeof (unsigned int) * ((cycleG->nodes->size) / 2));

    /* set all ddom_markers to zero */
    memset (ddom_markers, 0, sizeof (int) * nddom);

    /* set the number of dominoes traversed to zero */
    n_marked_dominoes = 0;

    /* set starting node and pid */
    current_node = start_node;
    current_pid = (start_node->id) / 2;

    /* mark the current pid */
    marray[current_pid] = 1;

    /* set odd-edge counter for current node */
    odd_num[current_pid] = n_marked_dominoes;

    /* choose next edge */
    if (!start_edge)
      current_edge = EGddpHeuristicSample (current_node,
                                           marray,
                                           n_marked_dominoes,
                                           odd_num, ddom_markers, 1);
    else
      current_edge = start_edge;

    do
    {

      /* if there are no edges left to choose from, we abort this path */
      if (!current_edge)
        break;

      /* update the current position and data structures */
      current_node = current_edge->head;
      current_pid = (current_node->id) / 2;
      node_data = (EGdijkstraNodeData_t *) (current_node->data);
      edge_data = (EGcycleData_t *) (current_edge->data);

      /* update dijkstra-node data */
      node_data->father = current_edge;

      /* update total path cost (from start_node) */
      path_val = EGdijkstraCostAdd (path_val, edge_data->weight);

      /* check if the current edge is a domino (odd) edge */
      if (edge_data->ddom)
        n_marked_dominoes++;

      /* forbid the edge is necessary */
      if (edge_data->ddom && ddom_markers)
        ddom_markers[edge_data->ddom->id] = 1;

      /* check if we have detected a cycle */
      if (marray[current_pid])
      {
        /* add the constraint if its not there already */
        /* note: this function checks if its violated  */
        rval =
          EGddpAddCutIfNew (ddpc_array, ddpc_dist, ddpc_angle_norm, ddpc_angle,
                            ddpc_size, ddpc_max_size, ddpc_best_angle,
                            ddpc_worst_angle, ddpc_worst_angle_id,
                            current_node, cycleG, dij_sol, marray, sarray, os,
                            1, mem);
        CHECKRVAL (rval);
        break;
      }
      /* check that we have not exceeded the bound */
      if ((ubound - EGdijkstraCostToLf (path_val)) < DDP_MIN_VIOLATION)
        break;
      /* mark the current pid */
      marray[current_pid] = 1;
      /* set odd-edge counter for current node */
      odd_num[current_pid] = n_marked_dominoes;
      /* choose next edge */
      current_edge = EGddpHeuristicSample (current_node, marray,
                                           n_marked_dominoes, odd_num,
                                           ddom_markers, 1);

    } while (1);
  }
  return 0;
}



int EGddpcComputeAll (EGmemPool_t * mem,
                      EGdGraph_t * bdG,
                      EGlist_t * dlist,
                      EGlist_t * ddpc_list,
                      int const nedges,
                      int k)
{
  register unsigned int i;
  int rval;
  size_t os[6];
  EGlistNode_t *n_it;
  EGdGraph_t *cycleG;
  EGdGraphNode_t *nleft,
   *nrigh;
  EGheap_t *my_heap;

  char *sarray;
  EGdijkstraCost_t max_val = EGdijkstraToCost (1 - DP_MAXVIOLATED_EPSILON);
  unsigned int *marray,
   *odd_num;
  double *max_weight;
  EGdGraphEdge_t **dij_sol;

  EGddpConstraint_t **ddpc_array;
  double *ddpc_dist;
  double *ddpc_angle_norm;
  double **ddpc_angle;
  unsigned int ddpc_size;
  unsigned int ddpc_max_size;
  double ddpc_best_angle;
  double ddpc_worst_angle;
  unsigned int ddpc_worst_angle_id;
  TEST (k != 1, "k (%d) is not one", k);

#if DDP_HEURISTIC
  EGtimer_t heur_timer;
  double max_ntime = 0.0;
#endif

#if DISPLAY_MODE
  fprintf (stdout, "DDPC: Computing dual DP constraints.\n");
#if DDP_HEURISTIC
  fprintf (stdout, "DDPC: Random Walk Heuristic: ON. Time Limit = %lf.\n",
           DDP_HEURISTIC_MAXTIME);
#endif
  fflush (stdout);
#endif

  os[EG_DIJ_DIST] = offsetof (EGdijkstraNodeData_t, dist);
  os[EG_DIJ_NDIST] = offsetof (EGdijkstraNodeData_t, ndist);
  os[EG_DIJ_FATHER] = offsetof (EGdijkstraNodeData_t, father);
  os[EG_DIJ_MARKER] = offsetof (EGdijkstraNodeData_t, marker);
  os[EG_DIJ_HCONNECTOR] = offsetof (EGdijkstraNodeData_t, hc);
  os[EG_DIJ_ELENGTH] = offsetof (EGcycleData_t, weight);

  cycleG = EGnewCycleGraph (mem, dlist, bdG, max_val);

  fprintf (stdout, "DDPC: Cycle graph has %d nodes and %d edges.\n",
           cycleG->nodes->size, cycleG->nedges);

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

  dij_sol = EGmemPoolSMalloc (mem, EGdGraphEdge_t *, cycleG->nedges);
  marray = EGmemPoolSMalloc (mem, unsigned int,
                             bdG->nodes->size);
  sarray = EGmemPoolSMalloc (mem, char,
                             bdG->nodes->size);
  odd_num = EGmemPoolSMalloc (mem, unsigned int,
                              bdG->nodes->size);
  max_weight = EGmemPoolSMalloc (mem, double,
                                 cycleG->nodes->size);

#if DDP_HEURISTIC

  /* set max_weight vector */
  rval = EGddpSetMaxWeight (cycleG, max_weight);
  CHECKRVAL (rval);

  /* set max-time per node */
  max_ntime = (DDP_HEURISTIC_MAXTIME / (cycleG->nodes->size / 2));

  /* re-set the heuristic timer */
  EGtimerReset (&heur_timer);

#endif

  ddpc_max_size = DDP_HEURISTIC_MAXCUTS;
  ddpc_array = EGmemPoolSMalloc (mem, EGddpConstraint_t *, ddpc_max_size);
  ddpc_dist = EGmemPoolSMalloc (mem, double,
                                ddpc_max_size);
  ddpc_angle_norm = EGmemPoolSMalloc (mem, double,
                                      ddpc_max_size);
  ddpc_angle = EGmemPoolSMalloc (mem, double *,
                                 ddpc_max_size + 1);
  ddpc_angle[ddpc_max_size] = EGmemPoolSMalloc (mem, double,
                                                ddpc_max_size);
  for (i = ddpc_max_size; i--;)
  {
    ddpc_angle[i] = EGmemPoolSMalloc (mem, double,
                                      ddpc_max_size);
    ddpc_angle[i][i] = 1.0;
  }
  ddpc_size = 0;


  for (n_it = cycleG->nodes->begin; n_it; n_it = ((n_it->next)->next))
  {

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

    nleft = (EGdGraphNode_t *) (n_it->this);
    nrigh = (EGdGraphNode_t *) (n_it->next->this);

    /* default: bad (exceeds the bound) */
    sarray[(nleft->id) / 2] = 'B';

    rval =
      EGpartialDijkstra (nleft, nrigh, EGdijkstraToCost (((double) 1.0)), os,
                         my_heap, cycleG);
    CHECKRVAL (rval);

    if (EGdijkstraGetMarker (nrigh, os) == UINT_MAX)
      continue;

    if ((1.0 - EGdijkstraCostToLf (EGdijkstraGetDist (nrigh, os))) <
        (DDP_MIN_VIOLATION))
      continue;

    rval =
      EGddpAddCutIfNew (ddpc_array, ddpc_dist, ddpc_angle_norm, ddpc_angle,
                        &ddpc_size, ddpc_max_size, &ddpc_best_angle,
                        &ddpc_worst_angle, &ddpc_worst_angle_id, nrigh,
                        cycleG, dij_sol, marray, sarray, os, 1, mem);
    CHECKRVAL (rval);

#if DDP_HEURISTIC
    EGtimerStart (&heur_timer);
    rval = EGddpAddHeuristicCuts_3 (ddpc_array, ddpc_dist, ddpc_angle_norm,
                                    ddpc_angle, &ddpc_size, ddpc_max_size,
                                    &ddpc_best_angle, &ddpc_worst_angle,
                                    &ddpc_worst_angle_id, nedges, nleft, 0, 1,
                                    1.0, cycleG, dij_sol, marray, sarray,
                                    odd_num, 0, 0, os, max_ntime, mem);
    CHECKRVAL (rval);
    EGtimerStop (&heur_timer);

#endif

  }

  /* pass the constraints from the heap to the list */
  fprintf (stderr, "***Heuristic Iterations: %llu\n", global_cnt);

  for (i = ddpc_size; i--;)
    EGlistPushBack (ddpc_list, ddpc_array[i]);

#if DISPLAY_MODE
  {
    EGddpConstraint_t *ddpc;
    double abs_slack,
      thresh,
      thresh_n;
    unsigned int range[DISPLAY_RESOLUTION],
      error = 0,
      maximally_violated = 0;
    double avg_teeth = 0.0;

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

    for (n_it = ddpc_list->begin; n_it; n_it = n_it->next)
    {
      ddpc = (EGddpConstraint_t *) (n_it->this);
      avg_teeth += ddpc->ndom;
      for (i = 0; i < DISPLAY_RESOLUTION; i++)
      {
        abs_slack = -1 * ddpc->slack_dp;
        thresh = 1 * (IDR) * i + 1 * (IDR);
        thresh_n = 1 * IDR * i;
        if (thresh_n <= abs_slack && abs_slack <= thresh)
        {
          range[i] += 1;
          break;
        }
      }

      if (i == DISPLAY_RESOLUTION)
        error += 1;

      if (-1 * ddpc->slack_dp > (1.0) - DDPCONSTRAINT_MAXVIOL_EPSILON)
        maximally_violated += 1;

    }
    avg_teeth = avg_teeth / ((double) ddpc_list->size);
#if DDP_HEURISTIC
    fprintf (stdout, "DDPC: Total heuristic running time: %lf seconds.\n",
             heur_timer.time);
#endif
    fprintf (stdout,
             "DDPC: Found %u dual DP-Constraints (avg number of teeth = %lf).\n",
             ddpc_list->size, avg_teeth);
    fprintf (stdout, "DDPC:    %u / %u are maximally violated.\n",
             maximally_violated, ddpc_list->size);
    fflush (stdout);

    fprintf (stdout, "DDPC: |slack values|\n");
    fflush (stdout);

    for (i = 0; i < DISPLAY_RESOLUTION; i++)
    {
      thresh = 1 * IDR * i + IDR * 1;
      thresh_n = 1 * IDR * i;
      if (range[i])
        fprintf (stdout, "[%8lf, %8lf]   %d\n", thresh_n, thresh, range[i]);
      fflush (stdout);
    }
    if (error)
      fprintf (stdout, "[%8lf,         )   %d    (out of range constraints)\n",
               1.0, error);
  }
#endif

#if DISPLAY_MODE
  fprintf (stdout, "DDPC: Finished finding dual DP cuts.\n");
  fflush (stdout);
#endif

  //EGfreeHeap (ddpc_heap, cycleG->mem);
  EGmemPoolSFree (ddpc_array, EGddpConstraint_t *, ddpc_max_size, mem);
  EGmemPoolSFree (ddpc_dist, double,
                  ddpc_max_size,
                  mem);
  EGmemPoolSFree (ddpc_angle_norm, double,
                  ddpc_max_size,
                  mem);
  for (i = ddpc_max_size + 1; i--;)
    EGmemPoolSFree (ddpc_angle[i], double,
                    ddpc_max_size,
                    mem);
  EGmemPoolSFree (ddpc_angle, double *,
                  ddpc_max_size + 1,
                  mem);

  EGmemPoolSFree (max_weight, double,
                  cycleG->nodes->size,
                  mem);
  EGmemPoolSFree (odd_num, unsigned int,
                  bdG->nodes->size,
                  mem);
  EGmemPoolSFree (marray, unsigned int,
                  bdG->nodes->size,
                  mem);
  EGmemPoolSFree (sarray, char,
                  bdG->nodes->size,
                  mem);
  EGmemPoolSFree (dij_sol, EGdGraphEdge_t *, cycleG->nedges, mem);

  EGfreeHeap (my_heap, mem);
  EGdGraphClearMP (cycleG, EGfreeCycleDataMP, EGfreeDijkstraNodeDataMP, 0, mem,
                   mem, 0);
  EGfreeDGraph (cycleG);


  return 0;
}

/* t = end node of the path, ddom_markers is work-array size of the number of
 * ddom. Checks if a ddom is used twice in the path. */
int EGddpIsSolutionValid (EGdGraphNode_t * t,
                          int *ddom_markers,
                          unsigned int num_ddom)
{

  EGdGraphEdge_t *e;
  EGdGraphNode_t *n;
  EGdijkstraNodeData_t *dndata;
  EGcycleData_t *cdata;

  unsigned int pos;

  memset (ddom_markers, 0, sizeof (int) * (num_ddom));

  n = t;
  pos = 0;

  while (((n->id) / 2 != (t->id) / 2) || !pos)
  {
    dndata = (EGdijkstraNodeData_t *) (n->data);
    e = dndata->father;
    cdata = (EGcycleData_t *) (e->data);

    if (cdata->ddom)
    {
      if (ddom_markers[cdata->ddom->id])
        return 0;
      else
        ddom_markers[cdata->ddom->id] += 1;
    }

    n = e->tail;

    pos += 1;
  }

  return 1;

}


int EGddpSetMaxWeight (EGdGraph_t * cycleG,
                       double *max_weight)
{

  EGdGraphEdge_t *e;
  EGdGraphNode_t *n;
  EGlistNode_t *e_it,
   *n_it;
  EGcycleData_t *cd;
  double val;

  for (n_it = cycleG->nodes->begin; n_it; n_it = n_it->next)
  {

    n = (EGdGraphNode_t *) (n_it->this);
    val = 0.0;

    for (e_it = n->out_edges->begin; e_it; e_it = e_it->next)
    {

      e = (EGdGraphEdge_t *) (e_it->this);
      cd = (EGcycleData_t *) (e->data);

      val += cd->pweight;

    }

    max_weight[(n->id)] = val;

  }

  return 0;

}

void EGddpcDisplay (EGddpConstraint_t * ddpc,
                    FILE * file)
{

  unsigned int i;

  fprintf (stderr, "F     [%u]: ", ddpc->nF);
  for (i = 0; i < ddpc->nF; i++)
    EGmengerDisplayEdgeBasic (ddpc->F[i], file);
  fprintf (stderr, "\n");

  fprintf (stderr, "ddom  [%u]: ", ddpc->ndom);
  for (i = 0; i < ddpc->ndom; i++)
    fprintf (stderr, "%d (1)(%lf) ", ddpc->ddom[i]->id,
             EGdijkstraCostToLf (ddpc->ddom[i]->value));
  fprintf (stderr, "\n");

  return;

}

---