mcf_cost_scaling1.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 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;
 }


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


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,
                                            int 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, class d_array>
bool acyclic_shortest_paths(const graph_t& G, rcost_t eps,
                                          const cap_array& cap,
                                          const flow_array& flow,
                                          const cost_array& cost,
                                          const pot_array& pot,
                                          d_array& dist,
                                          succ_array& succ,
                                          int& count)
{
  node_queue<graph_t,succ_array> Q(G,succ);
  dist_array indeg(G);

  int n = G.number_of_nodes();

  count = 0;

  node v;
  forall_nodes(v,G)
  { 
    //dist[v] = MAXINT;
    dist[v] = pot[v];

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

    indeg[v] = deg;

    if (deg  == 0) 
    { dist[v] = 0;
      Q.append(v);
     }
   }


  while (!Q.empty())
  { node u = Q.pop();
    rcost_t du = dist[u];
    n--;
    edge e;
    forall_out_edges(e,u)
    { node v = target(G,e,u);
      rcost_t rc = rcost(e,u,v);
      if (flow[e] == cap[e]) continue;
      rcost_t d = du + rc + eps;
      if (rc < 0 ) 
      { if (d < dist[v]) dist[v] = d;
        if (--indeg[v] == 0) Q.append(v);
       }
      else
        if (d < dist[v]) count++;
     }

    forall_in_edges(e,u)
    { node v = source(G,e,u);
      rcost_t rc = -rcost(e,v,u);
      if (flow[e] == 0) continue;
      rcost_t d = du + rc + eps;
      if (rc < 0 ) 
      { if (d < dist[v]) dist[v] = d;
        if (--indeg[v] == 0) Q.append(v);
       }
      else
        if (d < dist[v]) count++;
     }
   }

   return n == 0;
}



template<class flow_array, class cap_array, class cost_array>
bool bellman_ford_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,
                                               int n = -1)
{ 
  float tt = used_time();
  cout << "bellman_ford_refine: " << flush;

  if (n < 0) n = G.number_of_nodes();

  //dist_array dist(G);
  node_array<rcost_t,graph_t> dist(G,0);

  int count = 0;

  if (!acyclic_shortest_paths(G,eps,cap,flow,cost,pot,dist,succ,count)) 
  { cout << "not acyclic" << endl;
    return false;
  }


//  cout << string("#edges: %d",count) << endl;
//
//  if (count > G.number_of_edges()/2) 
//  { cout << string("too many edges: %d",count) << endl;
//    return false;
//   }

  node_queue<graph_t,succ_array> Q(succ);

  node v;
  forall_nodes(v,G) Q.append(v);

  node phase_last = Q.last();
  int  phase_count = 0;

  while (phase_count < n && !Q.empty())
  { 
    node u = Q.pop();
    rcost_t du = dist[u];

    edge e;
    forall_out_edges(e,u) 
    { if (flow[e] == cap[e]) continue;
      node v = target(G,e,u);    
      rcost_t d = du + rcost(e,u,v) + eps;
      if (d < dist[v])
      { dist[v] = d;
        if (!Q.member(v)) Q.append(v);
       }
     } 

    forall_in_edges(e,u) 
    { if (flow[e] == 0) continue;
      node v = source(G,e,u);    
      rcost_t d = du - rcost(e,v,u) + eps;
      if (d < dist[v])
      { dist[v] = d;
        if (!Q.member(v)) Q.append(v);
       }
     } 

    if (u == phase_last && !Q.empty()) 
    { phase_count++;
      phase_last = Q.last();
     }
  }


  cout << string("phase_count = %d  %.2f sec",phase_count,used_time(tt)) << endl;
  
  bool no_neg_cycle = false;

  if (Q.empty()) // no negative cycle
  { 
    // check
    forall_nodes(v,G)
    { edge e;
      forall_out_edges(e,v)
      { node w = target(G,e,v);
        if (flow[e] < cap[e]) 
          assert(dist[w] <= dist[v] + rcost(e,v,w) +eps); 
        if (flow[e] > 0) 
          assert(dist[v] <= dist[w] - rcost(e,v,w) +eps); 
       }
    }


    forall_nodes(v,G) 
    { rcost_t d = dist[v];
      pot[v] -= d;
     }

    no_neg_cycle = true;

    //check_eps_optimality(G,cost,pot,cap,flow,eps,"after bellman-ford");
  }

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


  return no_neg_cycle;
}




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

  forall_nodes(v,G) pot[v] += (dist[v]*epsd);
}