eg_menger.c

Go to the documentation of this file.

#include "eg_menger.h"


int EGmengerUpdateResidualGraph (EGdGraph_t * G,
                                 EGdGraphEdge_t * e,
                                 size_t * os)
{

  EGdGraphNode_t *swap;
  EGdijkstraCost_t cost;

  /* We unattach the edge from the graph.                                     */
  EGdGraphUnattachEdge (G, e);

  /* The cost changes sign.                                                   */
  cost =
    EGdijkstraCostMul (EGdijkstraToCost (-1.0), EGmengerGetEdgeCost (e, os));
  EGmengerSetEdgeCost (e, os, cost);

  /* We swap the head and tail.                                               */
  swap = e->head;
  e->head = e->tail;
  e->tail = swap;

  /* We re-attach the edge                                                    */
  EGdGraphAttachEdge (G, e);

  return 0;

}

int EGmengerRecoverGraphAndSolution (EGdGraph_t * G,
                                     size_t * os,
                                     EGdGraphNode_t * s,
                                     EGdGraphNode_t * t,
                                     unsigned int npath,
                                     EGdGraphEdge_t ** path,
                                     unsigned int *path_beg)
{

  unsigned int i,
    pos;
  EGlistNode_t *n_it,
   *e_it;
  EGdGraphNode_t *n;
  EGdGraphEdge_t *e;

  path_beg[0] = pos = 0;
  for (i = 0; i < npath; i++)
  {
    n = s;
    n_it = n->cn;
    while (n_it != t->cn)
    {
      for (e_it = n->in_edges->begin; e_it; e_it = e_it->next)
      {
        e = (EGdGraphEdge_t *) (e_it->this);
        if (EGmengerGetEdgeIsInSolution (e, os))
        {
          n = e->tail;
          n_it = n->cn;
          EGmengerSetEdgeIsInSolution (e, os, 0);
          EGmengerUpdateResidualGraph (G, e, os);
          path[pos++] = e;
          break;
        }
      }
      TEST (!e_it, "edge in solution not found for path %u", i);
    }
    path_beg[i + 1] = pos;
  }

  return 0;

}

/* EGmengerPathsADV: This algorithm is used to compute 'npaths' edge-disjoint
   paths between 's' and 't' having minimum total cost. The problem is tackled
   by solving a min-cost flow problem from 's' to 't' on a network with 
   edges of capacity one. 

   This program runs a succesive-shortest paths algorithm on G (as in Ahuja, 
   Magnanti, Orlin), taking advantage of the fact that for the problem we    
   need only consider capacity one edges.

   The output of the function is not the set of paths, however, but instead  
   is simply the final reduced graph and the value (menger_val) of the 
   optimal solution.                                                          */

/* This algorithm assumes that the fields 'isInSolution' of the edges are all 
   set to cero.                                                               */

int EGmengerPathsADV (EGdGraphNode_t * s,
                      EGdGraphNode_t * t,
                      EGdijkstraCost_t ubound,
                      unsigned int npaths,
                      unsigned int *nfpaths,
                      EGdijkstraCost_t * menger_val,
                      EGheap_t * dij_heap,
                      size_t * os,
                      EGdGraph_t * G)
{

  /* General use                                                              */
  int rval;
  unsigned int i;
  size_t dij_os[6];

  /* Algorithm-Variables                                                      */

  unsigned int f,
    path_nedges;
  EGdGraphNode_t *n;
  EGdGraphEdge_t *e;

  EGlistNode_t *n_it,
   *e_it;
  EGdijkstraCost_t prev_flow_val,
    dij_ubound,
    rcost,
    pi;

  /* Algorithm-Flags                                                          */
  *nfpaths = 0;

  /* Prepare the offsets which will be used by dijkstra                       */
  dij_os[EG_DIJ_DIST] = os[EG_MENGER_DIST];
  dij_os[EG_DIJ_NDIST] = os[EG_MENGER_NDIST];
  dij_os[EG_DIJ_FATHER] = os[EG_MENGER_FATHER];
  dij_os[EG_DIJ_MARKER] = os[EG_MENGER_MARK];
  dij_os[EG_DIJ_HCONNECTOR] = os[EG_MENGER_HEAP_CONNECTOR];
  dij_os[EG_DIJ_ELENGTH] = os[EG_MENGER_REDUCED_ECOST];

  /* Make sure that there are enough outbound edges in s.                     */
  if (s->out_edges->size < npaths)
    goto NO_SOLUTION;

  /* Make sure that there are enough inbound edges in t.                      */
  if (t->in_edges->size < npaths)
    goto NO_SOLUTION;

  /* Make sure the shortest distance from s to t isn't too large already.     */
  if (npaths * EGdijkstraCostToLf (EGmengerGetOrigDist (t, os)) >
      EGdijkstraCostToLf (ubound))
    goto NO_SOLUTION;

#if EG_MENGER_SAFEMODE > 1

/* do an initial check of negative-edges */

  for (n_it = G->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);
      TEST (EGdijkstraCostIsLess
            (EGmengerGetEdgeCost (e, os), EGdijkstraToCost (0.0)),
            "negative edge");
    }

  }

#endif

  /*                    ***** PART 1 *****
   * 
   * SEND THE FIRST UNIT OF FLOW FROM S TO T                       */

  /* Set the current value of the optimal solution.                           */
  *menger_val = prev_flow_val = EGmengerGetOrigDist (t, os);

  /* Set the dual prices on the node variables.                               */
  for (n_it = G->nodes->begin; n_it; n_it = n_it->next)
    EGmengerSetPi (n_it->this, os,
                   EGdijkstraCostMinus (EGmengerGetOrigDist (n_it->this, os)));

  /* Update the residual graph and mark the edges as 'taken'.                 */
  n = t;
  path_nedges = EGmengerGetOrigNdist (n, os);

  for (i = 0; i < path_nedges; i++)
  {
    e = EGmengerGetOrigFather (n, os);
    TEST (!e, "found node s (%u) before time %u / %u", n->id, i, path_nedges);
    n = e->tail;

    EGmengerUpdateResidualGraph (G, e, os);
    EGmengerSetEdgeIsInSolution (e, os, 1);
  }

  /* Update the reduced costs of the edges.                                   */
  for (n_it = G->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);
      /* rcost = cost[tail,head] + pi[head] - pi[tail]                                */
      rcost =
        EGdijkstraCostAdd (EGmengerGetEdgeCost (e, os),
                           EGmengerGetPi (e->head, os));
      rcost = EGdijkstraCostSub (rcost, EGmengerGetPi (e->tail, os));

#if EG_MENGER_SAFEMODE
      if (EGdijkstraCostToLf (rcost) < 0.0)
      {
        fprintf (stderr,
                 "negative edge, WARNING!   rcost = %lf  ecost = %lf  pi_h = %lf  pi_t = %lf\n",
                 EGdijkstraCostToLf (rcost),
                 EGdijkstraCostToLf (EGmengerGetEdgeCost (e, os)),
                 EGdijkstraCostToLf (EGmengerGetPi (e->head, os)),
                 EGdijkstraCostToLf (EGmengerGetPi (e->tail, os)));
        rcost = EGdijkstraToCost (0.0);
      }
#endif

      EGmengerSetEdgeReducedCost (e, os, rcost);
      TEST (EGdijkstraCostToLf (rcost) < 0.0,
            "\nERROR: Negative edge (%u,%u)! cost=%lf  pi_h = %lf  pi_t = %lf  val = %lf | %lf",
            e->head->id, e->tail->id,
            EGdijkstraCostToLf (EGmengerGetEdgeCost (e, os)),
            EGdijkstraCostToLf (EGmengerGetPi (e->head, os)),
            EGdijkstraCostToLf (EGmengerGetPi (e->tail, os)),
            EGdijkstraCostToLf (EGmengerGetEdgeReducedCost (e, os)),
            EGdijkstraCostToLf (rcost));
    }
  }

  /*                    ***** PART 2 *****
   * 
   * SEND THE OTHER UNITS OF FLOW FROM S TO T                      */

  /* We have found a path                                                     */
  *nfpaths = 1;
  //fprintf(stderr, "*** 0 val = %lf\n", EGdijkstraCostToLf(*menger_val));

  for (f = 1; f < npaths; f++)
  {

    /* We compute an upper bound, and run the partial dijkstra algorithm.     */

    dij_ubound = EGdijkstraCostDiv (EGdijkstraCostSub (ubound, *menger_val),
                                    EGdijkstraToCost ((double) (npaths - f)));
    dij_ubound = EGdijkstraCostSub (dij_ubound, prev_flow_val);

#if MENGER_USE_BELLFORD == 0
    rval = EGpartialDijkstra (s, t, dij_ubound, dij_os, dij_heap, G);
    CHECKRVAL (rval);
    if (EGdijkstraGetMarker (t, dij_os) == UINT_MAX)
      goto NO_SOLUTION;
#endif

#if MENGER_USE_BELLFORD
    rval = EGbellFord (s, dij_os, G);
    CHECKRVAL (rval);
#endif

    /* Check if we have marked 't' with a relevant value                      */

    /* Check if we have surpassed the upper bound.                            */
    if (EGdijkstraCostIsLess (dij_ubound, EGmengerGetDist (t, os)))
      goto NO_SOLUTION;

    /* Update the current value of the solution and the prev_flow variable.   */
    prev_flow_val = EGdijkstraCostAdd (prev_flow_val, EGmengerGetDist (t, os));
    *menger_val = EGdijkstraCostAdd (prev_flow_val, *menger_val);

    /* Update the dual variables: That is, for every marked node 'i' :
     * 
     * pi[i] = pi[i] + d[t] - d[i]                                            */

    for (n_it = G->nodes->begin; n_it; n_it = n_it->next)
      if (EGmengerGetMark (n_it->this, os) < UINT_MAX - 1)
      {
        pi =
          EGdijkstraCostAdd (EGmengerGetPi (n_it->this, os),
                             EGmengerGetDist (t, os));
        pi = EGdijkstraCostSub (pi, EGmengerGetDist (n_it->this, os));
        EGmengerSetPi (n_it->this, os, pi);
      }

    /* Update the residual graph and mark the edges as 'taken' or 'untaken'   */

    n = t;
    path_nedges = EGmengerGetNdist (n, os);
    for (i = 0; i < path_nedges; i++)
    {
      e = EGmengerGetFather (n, os);
      TEST (!e, "found node s (%u) before time %u / %u", n->id, i, path_nedges);
      n = e->tail;

      EGmengerSetEdgeIsInSolution (e, os,
                                   1 - EGmengerGetEdgeIsInSolution (e, os));
      EGmengerUpdateResidualGraph (G, e, os);
    }

    /* We have found another path                                             */
    *nfpaths += 1;

    /* Update the reduced cost of the edges                                   */
    for (n_it = G->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);
        rcost =
          EGdijkstraCostAdd (EGmengerGetEdgeCost (e, os),
                             EGmengerGetPi (e->head, os));
        rcost = EGdijkstraCostSub (rcost, EGmengerGetPi (e->tail, os));

#if EG_MENGER_SAFEMODE
        if (EGdijkstraCostToLf (rcost) < 0.0)
        {
          fprintf (stdout,
                   "WARNING! Negative reduced-cost edge (value = %lf).\n",
                   EGdijkstraCostToLf (rcost));
          rcost = EGdijkstraToCost (0.0);
        }
#endif

        EGmengerSetEdgeReducedCost (e, os, rcost);
        TEST (EGdijkstraCostToLf (rcost) < 0.0,
              "ERROR: Negative edge (%u,%u)! cost=%lf  pi_h = %lf  pi_t = %lf  val = %lf | %lf",
              e->head->id, e->tail->id,
              EGdijkstraCostToLf (EGmengerGetEdgeCost (e, os)),
              EGdijkstraCostToLf (EGmengerGetPi (e->head, os)),
              EGdijkstraCostToLf (EGmengerGetPi (e->tail, os)),
              EGdijkstraCostToLf (EGmengerGetEdgeReducedCost (e, os)),
              EGdijkstraCostToLf (rcost));
      }
    }

    //fprintf(stderr, "*** %u val = %lf\n", f, EGdijkstraCostToLf(*menger_val));

  }

  return 0;

NO_SOLUTION:

  *menger_val = EG_DIJKSTRA_COST_MAX;

  return 0;

}

int EGmengerEmergencyRecovery (EGdGraph_t * G,
                               size_t * os)
{

  unsigned int cnt = 0;

  EGlistNode_t *n_it,
   *e_it;

  EGdGraphNode_t *n;
  EGdGraphEdge_t *e;

  /* first, "unflip" solution edges */

FIX_GRAPH:

  cnt += 1;

  if (cnt > 2 * (G->nedges))
    return 1;

  for (n_it = G->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);
      if (EGmengerGetEdgeIsInSolution (e, os))
      {
        EGmengerSetEdgeIsInSolution (e, os, 0);
        EGmengerUpdateResidualGraph (G, e, os);
        goto FIX_GRAPH;
      }
    }
  }

  /* next, make sure all edges have non-negative weights */

  for (n_it = G->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);
      if (EGdijkstraCostToLf (EGmengerGetEdgeCost (e, os)) < 0.0)
        return 1;
    }
  }

  return 0;

}

---