mcf_cost_scaling_l.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 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 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,
                                          int& count)
{
  int n = G.number_of_nodes();

  b_queue<node> Q(n);

  dist_array indeg(G);

  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, int N=-1)
{ 
  float tt = used_time();
  cout << "bellman_ford_refine: " << flush;

  //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,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;
   }
*/

  int n = G.number_of_nodes();
  b_queue<node> Q(n+1);

  node_array<bool,graph_t> in_Q(G);


  int phase_count = 0;

  if (N < 0) N = n;

  node v;
  forall_nodes(v,G) 
  //forall(v,L) 
  { Q.append(v); 
    in_Q[v] = true;
   }


  Q.append((node)nil); // end marker

  while (phase_count < N && Q.size() > 1)
  { node u = Q.pop();

    if ( u == nil) 
    { phase_count++;
      Q.append((node) nil);
      continue;
     }
    rcost_t du = dist[u];

    in_Q[u] = false;

    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])
      { if (!in_Q[v])  Q.append(v);
        dist[v] = d;
        in_Q[v] = true;
      }
     } 

    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])
      { if (!in_Q[v])  Q.append(v);
        dist[v] = d;
        in_Q[v] = true;
      }
     } 

  }

  bool no_neg_cycle = false;
  
  if (Q.size() <= 1) // 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");
  }

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

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

template <class cap_array, class flow_array, class cost_array>
int adm_number(const graph_t& G, node v, 
                                    const cap_array& cap,
                                    const flow_array& flow,
                                    const cost_array& cost,
                                    const pot_array& pot)
{
   int count = 0;

   edge e;
   forall_out_edges(e,v)
   { node w = target(G,e,v);
     NT x = cap[e]-flow[e];
     if (x > 0)
     { rcost_t rc = rcost(e,v,w);
       if (rc < 0) count++;
     }
    }
   forall_in_edges(e,v)
   { node w = source(G,e,v);
     NT x = flow[e];
     if (x > 0)
     { rcost_t rc = -rcost(e,w,v);
       if (rc < 0) count++;
     }
    }

   return count;
}


template <class cap_array, class flow_array, class cost_array>
int adm_number(const graph_t& G, node v, 
                                    const cap_array& cap,
                                    const flow_array& flow,
                                    const cost_array& cost,
                                    const pot_array& pot,
                                    edge& min_e,
                                    node& min_s,
                                    rcost_t& min_rc,
                                    NT& min_x)
{
   min_rc = MAXINT;

   int count = 0;