mcf_cost_scaling4.t

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

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



template <class cap_array, class flow_array, class cost_array>
NT adm_cap(const graph_t& G, node v, const cap_array& cap,
                                        const flow_array& flow,
                                        const cost_array& cost,
                                        const pot_array& pot)
{
   NT adm_c = 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) adm_c += x;
     }
    }
   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) adm_c += x;
     }
    }

   return adm_c;
}



template <class cap_array, class cost_array, class flow_array>
bool cost_scaling(const graph_t& G, const cap_array& cap,
                                    const cost_array& cost,
                                    excess_array& excess,
                                    flow_array& flow)
{
  // returns true iff there exists a feasible flow

  int n = G.number_of_nodes();

  int global_update_delta     = (f_update > 0) ? int(f_update*n) : MAXINT;
  int global_update_threshold = global_update_delta;


  pot_array  pot;   // potential
  pot.use_node_data(G,0);

  succ_array succ;
  succ.use_node_data(G,0);

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

  NT C = 0;   // maximal cost of any edge in G

  cost_array& ca = (cost_array&)cost;

  edge e;
  forall_edges(e,G) 
  { flow[e] = 0;
    NT nc = n*cost[e];
    ca[e] = nc;
    if (nc > C) C = nc;
    else if (-nc > C) C = -nc;
   }


  int eps=1;
  while (eps < C) eps *= alpha;

  //int refinement_bound = eps/alpha;
  int refinement_bound = MAXINT;

  feasible_flow<NT,graph_t,succ_array,dist_array> ff;  

  if (!ff.run(G,excess,cap,flow)) return false;

/*
  if (beta > 1)
  { while (eps > 1 &&
            price_refinement(G,eps/beta,cap,flow,cost,pot,succ,eps)) 
      eps /= beta;
  }
*/

  //dist_array adm_array(G);

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

  counter push_count = 0;
  counter discharge_count = 0;

  int relabel_count = 0;

  while (eps > 1)
  {
   //float t1 = used_time();

   node v;

   // compute 0-optimal pseudoflow 

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

  double total_excess = 0;
  int    s_count = 0;

  forall_nodes(v,G) 
  { NT ev = excess[v];
    if (ev > 0) 
    { active.append(v);
      total_excess += ev;
      s_count++;
     }
    //adm_array[v] = 0;
   }

/*
   cout << string("eps = %9d  excess = %12.0f   |active| = %5d", 
                       eps, total_excess,s_count) << flush;
*/
    
   while (!active.empty())
   { 
     discharge_count++;

     node v = active.top();

     NT ev = excess[v];

     while (ev > 0)
     {
       rcost_t min_rc = MAXINT;
       edge min_e = 0;
       NT   min_x = 0;
  
       edge e;
       forall_out_edges(e,v)
       { NT f = flow[e];
         NT x = cap[e] - f;
         if (x == 0) continue; 
         node w = target(G,e,v);
         rcost_t rc = rcost(e,v,w);
         if (rc < 0)   // admissible
         { NT ew = excess[w];

#if defined(PUSH_LOOK_AHEAD)
           int adm_w = adm_array[w];
           if (adm_w == -1) adm_w = adm_cap(G,w,cap,flow,cost,pot);
           if (ew >= 0 && adm_w == 0)
           { pot[w] += eps;
             rc += eps;
             if (rc < min_rc) { min_rc = rc; min_e = e;  min_x = x; }
             if (f > ew) f = ew;
             if (f > 0) // push back
             { excess[w] -= f;
               flow[e] -= f;
               ev += f;
              }
              adm_array[w] = -1;
            }
           else
#endif
           { push_count++;
             if (ev <= x) 
             { flow[e] = f + ev;
               excess[w] = ew + ev;
               if (ew+ev > 0 && !active.member(w)) active.fast_append(w);
               ev = 0;
               goto FINISH_DISCHARGE;
              }
             flow[e] = f + x;
             excess[w] = ew + x;
             ev -= x;
             if (ew+x > 0 && !active.member(w)) active.fast_append(w);
            }
          }
         else
         if (rc < min_rc) { min_rc = rc; min_e = e; min_x = x; }
       }           
  
  
       forall_in_edges(e,v)
       { NT f = flow[e];
         if (f == 0) continue; 
         node w = source(G,e,v);
         rcost_t rc = -rcost(e,w,v);
         if (rc < 0) // admissible
         { NT ew = excess[w];

#if defined(PUSH_LOOK_AHEAD)
           int adm_w = adm_array[w];
           if (adm_w == -1) adm_w = adm_cap(G,w,cap,flow,cost,pot);
           if (ew >= 0 && adm_w == 0)
           { pot[w] += eps;
             rc += eps;
             if (rc < min_rc) { min_rc = rc; min_e = e; min_x = -f; }
             NT x = cap[e] - f;
             if (x > ew) x = ew;
             if (x > 0) // push back
             { excess[w] -= x;
               flow[e] += x;
               ev += x;
              }
              adm_array[w] = -1;
            }
           else
#endif
           { push_count++;
             if (ev <= f)
             { flow[e] = f - ev;
               excess[w] = ew + ev;
               if (ew+ev > 0 && !active.member(w)) active.fast_append(w);
               ev = 0;
               goto FINISH_DISCHARGE;
              }
             flow[e] = 0;
             excess[w] = ew + f;
             ev -= f;
             if (ew+f > 0 && !active.member(w)) active.fast_append(w);
            }
          }
         else
         if (rc < min_rc) { min_rc = rc; min_e = e; min_x = -f; }
       }
  
/*
       if (relabel_count > global_update_threshold && !active.empty())
       { global_price_update(G, cap, flow, excess, cost, pot, eps);
         global_update_threshold += global_update_delta;
         node u;
         forall_nodes(u,G) adm_array[u] = -1;
         continue;
        }
*/

       relabel_count++;
       pot[v] += (eps+min_rc);
  
       if (ev <= min_x || ev <= -min_x)
       { push_count++;
         if (min_x > 0)
         { node w = target(G,min_e,v);
           NT ew = excess[w];
           excess[w] = ew+ev;
           flow[min_e] += ev;
           if (ew+ev > 0 && !active.member(w)) active.fast_append(w);
          }
         else
         { node w = source(G,min_e,v);
           NT ew = excess[w];
           excess[w] = ew+ev;
           flow[min_e] -= ev;
           if (ew+ev > 0 && !active.member(w)) active.fast_append(w);
          }
         ev = 0;
        }

  
FINISH_DISCHARGE:;
  
     } // while ev > 0

     active.del_top(v);
     excess[v] = ev;
     //adm_array[v] = -1;

     if (relabel_count > global_update_threshold && !active.empty())
     { global_price_update(G, cap, flow, excess, cost, pot, eps);
       global_update_threshold += global_update_delta;
/*
       node u;
       forall_nodes(u,G) adm_array[u] = -1;
*/
      }

    } // while active not empty


   if (beta > 1 && eps < refinement_bound)
   {  while (eps/beta > 0 &&
             price_refinement(G,eps/beta,cap,flow,cost,pot,succ,eps)) 
        eps /= beta;
    }


   //check_eps_optimality(G,cost,pot,cap,flow,eps,"after eps-phase");

   eps /= alpha;

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

  } // while eps > 1                                    


  num_discharges = discharge_count;
  num_pushes = push_count;
  num_relabels = relabel_count;

  forall_edges(e,G) ca[e] /= n;

  return true;
}



public:

mcf_cost_scaling() : alpha(0), beta(0), f_update(0),
                     num_nodes(0), num_edges(0), num_discharges(0), 
                     num_pushes(0), num_relabels(0), num_updates(0), T(0) {}


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,
                           float f = 1.25,
                           int a = 18,
                           int b = 0)
{
  T = used_time();

  num_discharges = 0;
  num_relabels = 0;
  num_updates  = 0;

  alpha = a;
  beta  = b;
  f_update = f;

  num_nodes = G.number_of_nodes();
  num_edges = G.number_of_edges();


  excess_array excess;         // supply values after elimination of
  excess.use_node_data(G,0);   // lower capacity bounds

  ucap_array& cap = (ucap_array&)ucap;  // non-const reference to cap array
                                        // used to modify capacities
                                        // temporarily for elimination of
                                        // lower capacity bounds
  int lc_count = 0;

  node v;
  forall_nodes(v,G) excess[v] = supply[v];

  forall_nodes(v,G) 
  { edge e;
    forall_out_edges(e,v) 
    { NT lc = lcap[e];
      if (lc != 0)                  // nonzero lower capacity bound
      { node w  = target(G,e,v);
        excess[v] -= lc;
        excess[w] += lc;
        cap[e] -= lc;
        lc_count++;
       }
     }
  }
  
 
  bool feasible = cost_scaling(G, cap, cost, excess, flow);

  
  if (lc_count)
  { // adjust flow and restore capacities
    edge e;
    forall_edges(e, G) 
    { NT lc = lcap[e];
      if (lc != 0)
      { flow[e] += lc;
        cap[e] += lc;
       }
     }
   }

  T = used_time(T);

  if (feasible)
  { // check flow
    edge e;
    forall_edges(e,G)
    {
      if (flow[e] < 0 || flow[e] > cap[e])
           error_handler(1,"check_flow: illegal flow value");
     }
  
    forall_nodes(v,G)
    { NT ev = supply[v];
      edge e;
      forall_out_edges(e,v) ev -= flow[e];
      forall_in_edges(e,v) ev += flow[e];
      if (ev != 0)
           error_handler(1,"check_flow: non-zero excess");
     }
  }
    
  return feasible;
}


float cpu_time() { return T; }
 
void statistics(ostream& out)
{
  out << string("%5d nodes  %5d edges   a = %d  b = %d  f = %.2f",
         num_nodes,num_edges,alpha,beta,f_update) << endl;
  out << string("%8d discharges %10d pushes %8d relabels %3d updates",
                   num_discharges, num_pushes,num_relabels, num_updates) << endl;

  out << endl;
}


};


LEDA_END_NAMESPACE

#if LEDA_ROOT_INCL_ID == 450457
#undef LEDA_ROOT_INCL_ID
#include <LEDA/POSTAMBLE.h>
#endif