mcf_cost_scaling3.t

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

#if !defined(LEDA_ROOT_INCL_ID)
#define LEDA_ROOT_INCL_ID 450457
#include <LEDA/PREAMBLE.h>
#endif

#include <LEDA/graph_alg.h>
#include <LEDA/b_queue.h>
#include <LEDA/node_pq22.h>
#include <LEDA/std/math.h>
#include <LEDA/assert.h>

#include <LEDA/templates/feasible_flow.t>


/*
#define PUSH_LOOK_AHEAD
*/


#define rcost(e,v,w) (cost[e] - pot[v] + pot[w])


LEDA_BEGIN_NAMESPACE

/*
class counter  {
public:
  counter(int) {}
  void operator++(int) {}
  operator int() { return 0; }
};
*/

typedef int counter;


template <class graph, class succ_array>
class node_queue {

typedef typename graph::node node;
typedef typename graph::edge edge;

const node _sentinel;

node  _first;
node  _last;

succ_array& NEXT;

public:

node_queue(succ_array& succ) : _sentinel(node(0xffffffff)), NEXT(succ)
{ _first = _last = _sentinel; }


node_queue(const graph& G, succ_array& succ) : _sentinel(node(0xffffffff)),
                                               NEXT(succ)
{ node v;
  forall_nodes(v,G) NEXT[v] = NULL; 
  _first = _last = _sentinel;;
 }

bool empty() const { return _first == _sentinel; }

bool member(node v) const { return NEXT[v] != NULL; }

void append(node v)
{ 
  if (empty())
    _first = v;
  else
    NEXT[_last] = v;

  NEXT[v] = _sentinel;
  _last = v;
 }

void fast_append(node v) // precondition: not empty
{ 
  NEXT[_last] = v;
  NEXT[v] = _sentinel;
  _last = v;
 }


void fast_push(node v) // precondition: not empty
{ 
  NEXT[v] = _first;
  _first = v;
 }

void push(node v) 
{ NEXT[v] = _first;
  _first = v;
  if (_last == _sentinel) _last = v;
 }


node succ(node v) { return NEXT[v]; }

node top()       { return _first; }
node first()     { return _first; }
node last()      { return _last; }
node sentinel()  { return _sentinel; }


void del_top() 
{ node v = _first;
  _first = NEXT[v];
  NEXT[v] = NULL;
}

void del_top(node v) 
{ _first = NEXT[v];
  NEXT[v] = NULL;
}

node pop() 
{ node v = top();
  del_top(v);
  return v;
}


};




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<double,graph_t>,
                    class succ_array   = node_array<node,graph_t> >


class mcf_cost_scaling 
{
  typedef typename graph_t::node node;
  typedef typename graph_t::edge edge;
  typedef typename pot_array::value_type rcost_t;

  int   alpha;
  int   beta;
  float f_update;
  
  int num_nodes;
  int num_edges;
  int num_discharges;
  int num_pushes;
  int num_relabels;
  int num_updates;

  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);
}



template<class flow_array,class cap_array,class cost_array>
bool check_eps_optimality(const graph_t& G, const cost_array& cost,
                                            const pot_array& pot,
                                            const cap_array&   cap,
                                            const flow_array& flow,
                                            rcost_t eps, 
                                            string msg="")
{
  int count = 0;
  node v;
  forall_nodes(v,G)
  { edge e;
    forall_out_edges(e,v)
    { node w = target(G,e,v);
      if (flow[e] < cap[e] &&  rcost(e,v,w) < -eps 
       || flow[e] > 0      && -rcost(e,v,w) < -eps) count++; 
     }
   }

   if (count && msg != "")
         error_handler(1,"check_eps_opt: " + msg);

   return count == 0;
}


template<class cap_array, class flow_array, class cost_array>
bool price_refine_phase1(const graph_t& G, rcost_t eps,
                                          const cap_array& cap,
                                          const flow_array& flow,
                                          const cost_array& cost,
                                          const pot_array& pot,
                                          dist_array& dist,
                                          succ_array& succ)
{

  // sort admissible network topologically

  node_queue<graph_t,succ_array> Q(G,succ);

  int n = G.number_of_nodes();

  node v;
  forall_nodes(v,G)
  { 
    int deg = 0;

    edge e;
    forall_out_edges(e,v)
    { node w = target(G,e,v);
      if (flow[e] > 0  && -rcost(e,v,w) < 0) deg++;
     }

    forall_in_edges(e,v)
    { node w = source(G,e,v);
      if (flow[e] < cap[e] && rcost(e,w,v) < 0) deg++;
     }

    dist[v] = deg; // indeg[v]

    if (deg == 0) Q.append(v);
   }

  for(node u = Q.first(); u != Q.sentinel(); u = Q.succ(u))
  { n--; 
    edge e;
    forall_out_edges(e,u)
    { if (flow[e] == cap[e]) continue;
      node v = target(G,e,u);
      rcost_t rc = rcost(e,u,v);
      if (rc >= 0) continue;
      if (--dist[v] == 0) Q.append(v);
     }

    forall_in_edges(e,u)
    { if (flow[e] == 0) continue;
      node v = source(G,e,u);
      rcost_t rc = -rcost(e,v,u);
      if (rc >= 0 ) continue; 
      if (--dist[v] == 0) Q.append(v);
     }
   }

  if (n > 0) return false;  // cycle detected


  for(node u = Q.first(); u != Q.sentinel(); u = Q.succ(u))
  { int du = dist[u];
    edge e;
    forall_out_edges(e,u)
    { if (flow[e] == cap[e]) continue;
      node v = target(G,e,u);
      rcost_t rc = rcost(e,u,v);
      if (rc >= 0) continue;
      int d = du + int((rc+0.5)/eps);
      if (d < dist[v]) dist[v] = d;
     }

    forall_in_edges(e,u)
    { if (flow[e] == 0) continue;
      node v = source(G,e,u);
      rcost_t rc = -rcost(e,v,u);
      if (rc >= 0 ) continue; 
      int d = du + int((rc+0.5)/eps);
      if (d < dist[v]) dist[v] = d;
     }
   }

  return true;
}


template<class flow_array, class cap_array, class cost_array>
bool price_refinement(const graph_t& G, rcost_t eps, 
                                               const cap_array& cap,
                                               const flow_array& flow,
                                               const cost_array& cost,
                                               pot_array& pot, 
                                               succ_array& succ, rcost_t old_eps)

{
  float tt = used_time();
  cout << string("price refinement(%8.0f) ",eps) << flush;

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

  dist_array dist(G);

  node_pq22<int,graph_t> PQ(n+m);

  bool result = false;
  int count;

  while (price_refine_phase1(G,eps,cap,flow,cost,pot,dist,succ)) 
  { 
    // admissible network is acyclic

    int count = 0;

    node v;
    forall_nodes(v,G) 
    { int dv = dist[v];
      assert(dv <= 0);
      if (dv < 0) 
      { count++;
        PQ.insert(v,dv);
       }
     }

    if (count == 0) // eps-optimal
    { result = true;
      break;
     }

    while (!PQ.empty())
    { node v = PQ.del_min();
      int dv = dist[v];
      if (dv > 0) continue;
      dist[v] = 1-dv;              // mark v as scanned (dist > 0)
      edge e;
      forall_out_edges(e,v)
      { if (cap[e] == flow[e]) continue;
        node w = target(G,e,v); 
        if (dist[w] > 0) continue; // already scanned
        rcost_t rc = rcost(e,v,w);
        int d = dv;
        if (rc >= 0) 
           //d += int(std::ceil(rc/eps));
           //d += (int((rc+0.5)/eps)+1);
           d += int(rc/eps+1);
        if (dist[w] > d) 
        { dist[w] = d;
          PQ.insert(w,d);
         }
      }
      
      forall_in_edges(e,v)
      { if (flow[e] == 0) continue;
        node w = source(G,e,v); 
        if (dist[w] > 0) continue; // already scanned
        rcost_t rc = -rcost(e,w,v);
        int d = dv;
        if (rc >= 0) 
           //d += int(std::ceil(rc/eps));
           //d += (int((rc+0.5)/eps)+1);
           d += int(rc/eps+1);
        if (dist[w] > d) 
        { dist[w] = d;
          PQ.insert(w,d);
         }
      }
   }

  // increase potential
  forall_nodes(v,G) {
   int d = dist[v];
   if (d > 0) d = 1-d; // unmark v
   pot[v] -= (d*eps);
  }
 }

 if (result)
    cout << " succeeded ";
 else
    cout << " failed    ";

 node v;
 forall_nodes(v,G) succ[v] = 0;

 cout << string(" %.2f sec",used_time(tt)) << endl;

 return result;
}



      
    
    
template<class flow_array, class cap_array, class cost_array>
void global_price_update(const graph_t& G, const cap_array& cap,
                                           const flow_array& flow,
                                           const excess_array& excess,
                                           const cost_array& cost,
                                           pot_array& pot,
                                           rcost_t epsd)
{

 num_updates++;

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

 int inf = alpha*n;

 node_pq22<int,graph_t> PQ(n+m);

 dist_array dist(G);

 node v;
 forall_nodes(v,G) 
 { NT ev = excess[v];
   if (ev < 0)
   { dist[v] = 0;
     PQ.insert(v,0);
    }
   else 
     dist[v] = inf;
  }

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

    edge e;
    forall_out_edges(e,u)
    { if (flow[e] == 0) continue;
      node v = target(G,e,u);
      int dv = (int)dist[v];
      if (dv <= du) continue;
      double rc = -rcost(e,u,v);
      int c = du;
      if (rc >=0) c += int(rc/epsd)+1;
      if (c < 0) continue; // overflow
      if (c < dv)  
      {  PQ.insert(v,c);
         dist[v] = c;
       }
     }

    forall_in_edges(e,u)
    { if (flow[e] == cap[e]) continue;
      node v = source(G,e,u);
      int dv = (int)dist[v];
      if (dv <= du) continue;
      double rc = rcost(e,v,u);
      int c = du;
      if (rc >=0) c += int(rc/epsd)+1;
      if (c < 0) continue; // overflow
      if (c < dv)  
      {  PQ.insert(v,c);
         dist[v] = c;
       }
     }
   }