t

A Library of Efficient Data Types and Algorithms,封装了常用的ADT及其相关算法的软件包

#include <LEDA/graph_alg.h>
#include <LEDA/node_pq22.h>
#include <LEDA/b_stack.h>
#include <LEDA/assert.h>

const unsigned long mask1 = 0x80000000;
const unsigned long mask2 = 0x7fffffff;

LEDA_BEGIN_NAMESPACE

template <class NT, class graph_t      = graph,
                    class pot_array    = node_array<NT,graph_t>,
                    class excess_array = node_array<NT,graph_t>,
                    class dist_array   = node_array<NT,graph_t>,
                    class pred_array   = node_array<edge,graph_t>,
                    class cap_array    = edge_array<NT,graph_t>,
                    class node_pq_t    = node_pq22<NT,graph_t> >

class mcf_cap_scaling 
{
  typedef typename graph_t::node node;
  typedef typename graph_t::edge edge;
  typedef node_pq_t node_pq;
  
  float T;

node source(const graph_t& G, edge e, node t)
{ if (graph_t::category() == opposite_graph_category)
    return G.opposite(e,t);
  else
    return G.source(e);
}
  
node target(const graph_t& G, edge e, node s)
{ if (graph_t::category() == opposite_graph_category)
    return G.opposite(e,s);
  else
    return G.target(e);
}



inline void flip_sign(const NT& x) { (NT&)x = -x; }



template<class cost_array, class flow_array>

node DIJKSTRA_IN_RESIDUAL(const graph_t& G, node s, int delta, 
                             const cost_array&   cost,
                             const cap_array&    cap,
                             const flow_array&   flow,
                             const excess_array& excess,
                             pot_array&  pi,
                             dist_array& dist,
                             pred_array& pred,
                             node_pq22<NT,graph_t>&  PQ)

{

  int n = G.number_of_nodes();
  b_stack<node> labeled(n);   // list of all (permanently  or temporarily)
                              // labeled nodes

  node t = nil;
  NT dt = MAXINT;

  dist[s] = 0;
  PQ.insert(s,0);


  while (!PQ.empty())
  { 
    node u = PQ.del_min(dist);
    NT  du = dist[u];

    if (du == -1) continue;
    
    labeled.push(u);

    du += pi[u]; 
    pi[u] = du;
    dist[u] = -1;

    if (u == t) break;

    edge e;
    forall_out_edges(e,u) 
      if (cap[e]-flow[e] >= delta)  // e in delta-residual graph
      { node v = target(G,e,u);
        NT c = du - pi[v] + cost[e];

        if (c >= dt) continue;
        if (excess[v] <= -delta) { t = v; dt = c; }

        if (c < dist[v])
        { PQ.insert(v,c); 
          dist[v] = c;
          pred[v] = e;
         }
       }

    forall_in_edges(e,u) 
      if (flow[e] >= delta)  // e in delta-residual graph
      { node v = source(G,e,u);
        NT c = du - pi[v] - cost[e];

        if (c >= dt) continue;
        if (excess[v] <= -delta) { t = v; dt = c; }

        if (c < dist[v])
        { PQ.insert(v,c); 
          dist[v] = c;
          pred[v] = edge(((unsigned long)e) | mask1);
         }
       }

   }


  while (!labeled.empty())
  { node v = labeled.pop();
    if (t) pi[v] -= dt;
    dist[v] = MAXINT;
   }

  for(int i=0; i<PQ.size(); i++) 
  { node v = PQ[i];
    dist[v] = MAXINT;
   }

  PQ.clear();

  return t;
}



public:


template <class lcap_array, class ucap_array, class cost_array, 
          class supply_array, class flow_array>
bool run(const graph_t& G, 
                           const lcap_array& lcap,
                           const ucap_array& ucap,
                           const cost_array& cost,
                           const supply_array& supply,
                           flow_array& flow, int f = 16)
{
  T = used_time();

  int n = G.number_of_nodes();
  int m = G.number_of_edges();

  int count = 0;		// number of augmentations
  node v;
  edge e;

  NT U = 0;		// maximal supply/demand or capacity

  forall_edges(e, G) { 
    assert(cost[e] >= 0);
    if (ucap[e] > U)   U = ucap[e];
  }

  pot_array  pi;
  pi.use_node_data(G,0);

  excess_array  excess;
  excess.use_node_data(G,0);

  dist_array dist;
  dist.use_node_data(G,MAXINT);

  pred_array pred;
  pred.use_node_data(G,0);

  cap_array cap;
  cap.use_edge_data(G,0);


  node_pq PQ(n+m);


  forall_nodes(v, G) 
  { NT s = supply[v];
    if ( s > U) U =  s; 
    if (-s > U) U = -s;
    excess[v] = s;
   }
  

  list<edge> neg_edges;		// list of negative cost edges
  int lc_count = 0;             // number of low-cap edges

  forall_nodes(v,G)
  { edge e;
    forall_out_edges(e,v)
    { flow[e] = 0;
      cap[e] = ucap[e];
      if (cost[e] < 0) neg_edges.append(e);

      // removing nonzero lower bounds (Ahuja/Magnanti/Orlin, section 2.4,
      // p. 39 [fig. 2.19 is not correct!]
        
      NT lc = lcap[e];
  
      if (lc != 0) 
      {	// nonzero lower capacity
        lc_count++;
        node i = v;
        node j = target(G,e,v);
        excess[i] -= lc;
        excess[j] += lc;
        cap[e] -= lc;
       }
     }
   }



  // eliminate negative cost edges by edge reversals 
  // (Ahuja/Magnanti/Orlin, section 2.4, p. 40)

/*
  forall(e, neg_edges) 
  { node u = G.source(e);
    node v = G.target(e);
    excess[u] -= cap[e];
    excess[v] += cap[e];
    flip_sign(cost[e]);
    G.rev_edge(e);
   }
*/


  int delta = 1;

  U = (U*m)/(n*n);

  while (U >= delta) delta *= 2;

  delta *= f;

  while (delta > 0) 
  { // delta scaling phase

    // Saturate all edges entering the delta residual network which have
    // a negative reduced edge cost. Then only the reverse edge (with positive
    // reduced edge cost) will stay in the residual network. 

  node u;
  forall_nodes(u,G)
  { edge e;
    forall_out_edges(e,u)
    { NT fe = flow[e];
      if (fe >= delta)
      { node v = target(G,e,u);
        NT c = cost[e] + pi[u] - pi[v];
        if (c > 0)
        { excess[u] += fe;
          excess[v] -= fe;
          flow[e] = 0;
          continue;
         }
       }
      NT ce  = cap[e];
      NT cmf = ce-fe;
      if (cmf >= delta)
      { node v = target(G,e,u);
        NT c = cost[e] + pi[u] - pi[v];
        if (c < 0)
        { excess[u] -= cmf;
          excess[v] += cmf;
          flow[e] = ce;
         }
       }
     }
   }


    // As long as there is a node s with excess >= delta compute a minimum
    // cost augmenting path from s to a sink node t with excess <= -delta in 
    // the delta-residual network and augment the flow by at least delta units 
    // along P (if there are nodes with excess > delta and excess < -delta
    // respectively

  node s;
  
  forall_nodes(s,G)
  { NT es = excess[s];

    while (es >= delta)
    { 
      count++;
  
      node t = DIJKSTRA_IN_RESIDUAL(G,s,delta,cost,cap,flow, excess,pi,
                                    dist,pred,PQ);

      if (t == nil) break; // there is no node with excess < -delta

      NT et = -excess[t];
  
      // compute maximal amount f of flow that can be pushed through P

      NT f = min(es,et);

      for(v = t; v != s; v = G.opposite(e,v))
      { e = pred[v];
        NT rc;
        if (((unsigned long)e) & mask1)
          { e = edge(((unsigned long)e) & mask2);
            rc = flow[e];
           }
        else
           rc = cap[e]-flow[e];

        if (rc < f) f = rc;
       }
  
      // add f units of flow along the augmenting path 

      for(v = t; v != s; v = G.opposite(e,v))
      { e = pred[v];
        if (((unsigned long)e) & mask1)
          { e = edge(((unsigned long)e) & mask2);
            flow[e] -= f;
           }
        else
           flow[e] += f;
       }
  
       es -= f;
       et -= f;

       excess[t] = -et;
  
    } // end of delta-phase

    excess[s] = es;
  }
  
  delta /= 2;
  
 } // end of scaling
  


/*
  // restore negative edges
  edge x;
  forall(x, neg_edges) 
  { flow[x] = cap[x] - flow[x];
    G.rev_edge(x);
    flip_sign(cost[e]);
   } 
*/


  // adjust flow

  if (lc_count > 0)
  { edge x;
    forall_edges(x,G) flow[x] += lcap[x];
   }

  bool feasible = true;

  forall_nodes(v,G)
    if (excess[v]  != 0) feasible = false;

  T = used_time(T);

  return feasible;
}


float cpu_time() const { return T; }

};



LEDA_END_NAMESPACE