mcf_cost_scaling_l.t

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

   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++;
       else if (rc < min_rc) { min_rc = rc; min_x = x; min_e = e; min_s = v; }
     }
    }
   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++;
       else if (rc < min_rc) { min_rc = rc; min_x = x; min_e = e; min_s = w; }
     }
    }

   return count;
}



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;

  int beta4 = beta*beta*beta*beta;

  pot_array  pot;   // potential
  pot.use_node_data(G,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;


  feasible_flow<NT,graph_t> ff;  
  if (!ff.run(G,excess,cap,flow)) return false;

  if (beta > 1)
  { while (bellman_ford_refinement(G,eps/beta,cap,flow,cost,pot,n/100)) 
    { eps /= beta;
      cout << "new eps0 = " << eps << endl;
    }
  }


  b_queue<node> active(n);

  counter push_count = 0;
  counter discharge_count = 0;
  counter relabel_count = 0;


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

   // compute 0-optimal pseudoflow 

   node v;
   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;

  forall_nodes(v,G) 
  { NT ev = excess[v];
    if (ev > 0) 
    { active.append(v);
      total_excess += ev;
     }
//assert(adm_number(G,v,cap,flow,cost,pot) == 0);
   }

/*
   cout << string("eps = %9d  excess = %12.0f   |active| = %5d", 
                       eps, total_excess,active.size()) << flush;
*/

  //node_array<int,graph_t> adm_count(G,0);

   pot[active.top()] += eps;

   while (!active.empty())
   { 
     discharge_count++;

     node v = active.top();


     NT ev = excess[v];

     rcost_t min_rc = MAXINT;
     edge min_e = 0;
     node min_s = 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)  
       { NT ew = excess[w];
         if (ew >= 0 && adm_number(G,w,cap,flow,cost,pot) == 0)
         { 
           pot[w] += eps;
           rc += eps;
           if (rc < min_rc) { min_rc = rc; min_e = e; min_s = v; min_x = x; }

           if (ew > 0 && f > 0)
           { // push back
             if (ew > f)
             { excess[w] -= f;
               flow[e] = 0;
               ev += f;
              }
            }
          }
         else
         { push_count++;
           if (ev <= x) 
           { flow[e] = f + ev;
             excess[w] = ew + ev;
             if (ew <= 0 && ev > -ew) active.append(w);
             goto NEXT_NODE;
            }
           flow[e] = f + x;
           excess[w] = ew + x;
           ev -= x;
           if (ew <= 0 && x > -ew) active.append(w);
          }
        }
       else
       if (rc < min_rc) { min_rc = rc; min_e = e; min_s = v; 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)
       { NT ew = excess[w];
         if (excess[w] >= 0 && adm_number(G,w,cap,flow,cost,pot) == 0)
         { pot[w] += eps;
           rc += eps;
           if (rc < min_rc) { min_rc = rc; min_e = e; min_s = w; min_x = f; }

           NT x = cap[e] - flow[e];
           if (ew > 0 && x > 0)
           { // push back
             if (ew > x)
             { excess[w] -= x;
               flow[e] += x;
               ev += x;
              }
            }
          }
         else
         { push_count++;
           if (ev <= f)
           { flow[e] = f - ev;
             excess[w] = ew + ev;
             if (ew <= 0 && ev > -ew) active.append(w);
             goto NEXT_NODE;
            }
           flow[e] = 0;
           excess[w] = ew + f;
           ev -= f;
           if (ew <= 0 && f > -ew) active.append(w);
          }
        }
       else
       if (rc < min_rc) { min_rc = rc; min_e = e; min_s = w; min_x = f; }
      }

     excess[v] = ev;


     if (relabel_count > global_update_threshold && !active.empty())
     { global_price_update(G, cap, flow, excess, cost, pot, eps);
       global_update_threshold += global_update_delta;
       continue;
      }

     relabel_count++;
     pot[v] += (eps+min_rc);

     if (ev <= min_x) 
     { push_count++;
       if (v == min_s)
       { node w = target(G,min_e,v);
         NT ew = excess[w];
         NT ew_new = ew + ev;
         if (ew <= 0 && ew_new > 0) active.append(w);
         excess[w] = ew_new;
         flow[min_e] += ev;
        }
       else
       { node w = source(G,min_e,v);
         NT ew = excess[w];
         NT ew_new = ew + ev;
         if (ew <= 0 && ew_new > 0) active.append(w);
         excess[w] = ew_new;
         flow[min_e] -= ev;
        }
       goto NEXT_NODE;
      }
    
     continue;

NEXT_NODE:;
     excess[v] = 0;
     active.del_top();


    } // while active not empty


   if (eps < beta4)
   { while (eps > 1 && 
            bellman_ford_refinement(G,eps/beta,cap,flow,cost,pot,n/100)) 
     { eps /= beta;
       cout << "new eps = " << eps << endl;
     }
   }

   //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