mincostflow.t

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

/*******************************************************************************
+
+  LEDA 4.5  
+
+
+  mincostflow.t
+
+
+  Copyright (c) 1995-2004
+  by Algorithmic Solutions Software GmbH
+  All rights reserved.
+ 
*******************************************************************************/

// $Revision: 1.5 $  $Date: 2004/02/06 11:20:20 $

#ifndef MINCOSTFLOW_H
#define MINCOSTFLOW_H

#include <LEDA/node_slist.h>

//#include <LEDA/std/limits.h>
//#include <LEDA/std/iostream.h>

LEDA_BEGIN_NAMESPACE


#ifdef  TIMER_ENABLED    
#define START_TIMER      float T = used_time();
#define STOP_TIMER(time) time += used_time(T);
#else
#define START_TIMER      
#define STOP_TIMER(time) 
#endif


// hack for static graphs

template <class node>
inline GenPtr First_Adj_Edge(node v, int i)
{ if (i == 0) return v->first_out_edge();
  return v->first_in_edge();
}


template <class NT, class graph = leda::graph,
                    class cap_array = edge_array<NT,graph>,
                    class exc_array = node_array<NT,graph>,
                    class pot_array = node_array<double,graph>,
                    class cur_array = node_array<typename graph::edge,graph>,
                    class act_list  = base_node_slist<graph,node_array<typename graph::node,graph> > >
class mincostflow
{     
  int num_discharges[16];
  int num_relabels[16];  
  int num_pushes[16];
  int num_refines;

  float epsilon;
  float scaling_factor;

  float init_time;  
  float refine_time;  
  float push_time;
  float relabel_time;
  float seek_time;
  float check_time;
  float discharge_time;  
  float total_time;
  
  NT     total_flow;
  double total_cost;
  
  int num_nodes;
  int num_edges;
  
  public:
    
  typedef typename graph::node node;
  typedef typename graph::edge edge;
 
  mincostflow(float factor = 8) : scaling_factor(factor) { reset_counter(); }
   
  template <class cost_array, class lcap_array, class ucap_array,
                              class sup_array,  class flow_array>
  void run(const graph& G, const cost_array& cost, const lcap_array& lcap, 
           const ucap_array& ucap, const sup_array& supply, flow_array& flow)
  {      
    float T = used_time();    
        
    cap_array cap;    
    exc_array exc; 
    cur_array cur;    

    init(G,cost,lcap,ucap,supply,cap,exc,cur); 
    
    pot_array pot;    

    flow.use_edge_data(G,0);
    pot.use_node_data(G,0);
        
    act_list active(G);
    
    double bound = 1 / double(num_nodes);  
    
    while (epsilon > bound)
    { refine(G,cost,cap,flow,exc,pot,cur,active);           
      reduce_eps();                       
    }     

    total_flow = 0;
    
    edge e;
    forall_edges(e,G)
    { flow[e]    += lcap[e];
      total_flow += flow[e];
      total_cost += flow[e] * cost[e];
    }

    total_time = used_time(T);
  }


  void statistics(ostream& os = cout)
  {     
    int sum[4], i;
    for (i = 0; i < 4; i++) sum[i] = 0;      
      
    for (i = 0; i < num_refines; i++)
    { sum[0] += num_discharges[i]; 
      sum[1] += num_relabels[i];
      sum[2] += num_pushes[i];  
    }
    
    os << "\n";
    os << "c cost scaling version 1.0\n";
    os << "c\n";
    os << string("c nodes: %15.0d     edges: %15.0d\n", num_nodes, num_edges);  
    os << string("c scaling-factor:  %5.0f\n", scaling_factor);
    os << "c\n";
    os << string("c time:  %15.2f     cost:  %15.0f\n", total_time, total_cost); 
    os << string("c refines:    %10.0d     discharges: %10d\n", num_refines, sum[0]);
    os << string("c pushes:     %10.0d     relabels:   %10d\n", sum[2], sum[1]);
    os << "c\n"; 
    os << "c discharge     pushes   relabels\n"; 
          
    for (i = 0; i < num_refines; i++)
      os << string("c%10d %10d %10d\n", num_discharges[i], num_pushes[i], num_relabels[i]);    
    
    os << "c\n";    
    
#ifdef TIMER_ENABLED    
    os << string("c init:      %10.2f\n", init_time);
    os << string("c refines:   %10.2f     discharges:        %10.2f\n", refine_time, discharge_time);
    os << string("c pushes:    %10.2f     relabels:          %10.2f\n", push_time, relabel_time);
    os << string("c cur. edge: %10.2f     admissible checks: %10.2f\n", seek_time, check_time);
    os << "\n";
#endif  
  }
  
  NT     flow()     const { return total_flow; }
  double cost()     const { return total_cost; }  
  float  cpu_time() const { return total_time; }

  private:   
  
  // ----------------------------------------------------------------------
  // check epsilon-optimality:
  // for each edge (v,w): cp(v,w) >= -epsilon && rcap(v,w) > 0
  // for each edge (v,w): cp(v,w) < epsilon -> rcap(v,w) == 0
  // ----------------------------------------------------------------------

  template <class cost_array, class flow_array>
  bool check(const graph& G, const cost_array& cost, const cap_array& cap, 
             const flow_array& flow, const pot_array& pot)
  {
    typedef typename pot_array::value_type PT;    
      
    edge e;
    forall_edges(e,G)
    {
      node v = G.source(e);
      node w = G.target(e);
      PT cp = cost[e] + pot[v] - pot[w];
      if (cp >= -epsilon)
      { if (flow[e] < cap[e] || cap[e] == 0) continue;
        cerr << "\nmincostflow::check_eps_optimality() \n";
        cerr << "flow is not epsilon-optimal:\n";
        cerr << "epsilon = " << epsilon << "\n";
        cerr << "cp(v,w) = " << cp << "\n";
        cerr << "rc(v,w) = " << cap[e] - flow[e] << "\n";
        return false;
      }
      else
      { if (flow[e] == cap[e]) continue;
        cerr << "\nmincostflow::check_eps_optimality() \n";
        cerr << "flow is not epsilon-optimal:\n";
        cerr << "epsilon = " << epsilon << "\n";
        cerr << "cp(v,w) = " << cp << "\n";
        cerr << "rc(v,w) = " << cap[e] - flow[e] << "\n";
        return false;
      }

      return true;
    }  
  }

  
  void reduce_eps() { epsilon /= scaling_factor; }

  // ----------------------------------------------------------------------
  // refine:
  // determines an eps. optimal flow by performing discharge operations
  // on all active nodes (nodes with excess > 0)
  // ----------------------------------------------------------------------

  template <class cost_array, class flow_array, class node_slist>
  void refine(const graph& G, const cost_array& cost, const cap_array& cap, 
                    flow_array& flow, exc_array& excess, pot_array& pot,
                    cur_array& cur, node_slist& active)
  { START_TIMER;
    
    // -----------------------------------------------------------------
    // converts the flow into an eps-optimal pseudoflow 
    // -----------------------------------------------------------------
      
    edge e;
    forall_edges(e,G)
    { node v = G.source(e);
      node w = G.target(e);
      
      if (cost[e] + pot[v] - pot[w] < -epsilon)
      { NT df = cap[e] - flow[e];
        excess[v] -= df;
        excess[w] += df;
        flow[e]   += df;          
      }
      else
      { NT df = flow[e];
        excess[v] += df;
        excess[w] -= df;
        flow[e]    = 0;
      }
    }
    
    // -----------------------------------------------------------------
    // insert nodes with excess > 0 into node queue active
    // -----------------------------------------------------------------

    node v;
    forall_nodes(v,G)
      if (excess[v] > 0) active.append(v);   
    
    while (!active.empty())
    { node v = active.pop();         
      if (excess[v] > 0) 
        discharge(G,cost,cap,flow,excess,pot,cur,active,v);
    }  
        
    num_refines++;
    STOP_TIMER(refine_time);
  }



  // ------------------------------------------------------------
  // push(v,w)
  // applicability: v is active, 
  //                residual capacity of e: rescap(e) > 0,
  //                reduced cost of e: cp(e) < 0 
  //                   
  // action: send min(exc[v],rescap(e)) units flow from v to w
  // ------------------------------------------------------------

  template <class flow_array, class node_slist>
  void push(const graph& G, const cap_array& cap, 
            flow_array& flow, exc_array& excess, 
            node_slist& active, node v, node w, edge e)
  { START_TIMER;

    NT df = excess[v];
    
    if (w == G.target(e)) // forward
    { if (cap[e] - flow[e] < df) df = cap[e] - flow[e];
      flow[e] += df;
    }
    else                  // backward
    { if (flow[e] < df) df = flow[e];
      flow[e] -= df;    
    }