mincostflow.t

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

    
    excess[v] -= df;
    excess[w] += df;
    
    if (!active.member(w)) active.append(w);          

    num_pushes[num_refines]++;          
    STOP_TIMER(push_time);
  }



  // ------------------------------------------------------------
  // push_back(w,v) used for push-lookahead heuristic
  // applicability: residual capacity of (w,v) > 0
  //                reduced cost of (w,v) < 0
  //
  // action: send min(exc[w],rescap(w,v)) units flow from w to v
  // -----------------------------------------------------------

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

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


  // ---------------------------------------------------------
  // is_admissible(current_edge(v))
  // checks whether the current edge of v is admissible;
  // two conditions must be accomplished:
  // 1. residual capacity of e: rescap(e) > 0
  // 2. reduced cost of e: cp(e) < 0
  // ---------------------------------------------------------

  template <class cost_array, class flow_array>
  bool is_admissible(const graph& G, const cost_array& cost, 
                     const cap_array&  cap, const flow_array& flow,
                     const pot_array&  pot, const cur_array& cur, node v)
  { START_TIMER;
    
    typedef typename pot_array::value_type PT;    
  
    NT rc;
    PT cp;  

    edge e = cur[v];
    node w = G.opposite(e,v);
    if (w == G.target(e))    // forward
    { rc = cap[e] - flow[e];
      cp = cost[e] + pot[v] - pot[w];
    }
    else   
    { rc = flow[e];          // backward
      cp = -cost[e] + pot[v] - pot[w];
    }

    STOP_TIMER(check_time);
    return rc > 0 && cp < 0; 
  }

  
  // -----------------------------------------------------------
  // relabel(v)
  // applicability: v is active and
  //                v hasn't an admissible edge
  //
  // action: replace pot(v) by max(pot(w) - cost(e) - epsilon)
  //         for all in- and out-edges e of v
  //
  //  (v,w) : pot[v] = pot[w] - cost(v,w) - epsilon
  //  (w,v) : pot[v] = pot[w] - (-cost(w,v)) - epsilon
  //
  // note: there is an admissible edge of v, we set e 
  //       to the current edge of v and leave the function
  //       without relabeling v
  // ----------------------------------------------------------

  template <class cost_array, class flow_array>
  bool relabel(const graph& G, const cost_array& cost, 
               const cap_array& cap, const flow_array& flow, 
               exc_array& excess, pot_array& pot, cur_array& cur, node v)
  { START_TIMER;

    typedef typename pot_array::value_type PT;

    PT min_p = -MAXINT; //-numeric_limits<PT>::max();
    PT cur_p = min_p;
    
    edge e, cur_e = 0;
    forall_out_edges(e,v)
    { if (flow[e] < cap[e])          // residual cap(e) > 0
      { node w = G.target(e);
        PT  dp = pot[w];
      
        if (cost[e] + pot[v] < dp)   // reduced cost cp(v,w) < 0
        { cur[v] = e;            
          STOP_TIMER(seek_time);
          return false;                   
        }
   
        dp -= cost[e];
   
        if (dp > cur_p) 
        { cur_p = dp;
          cur_e = e;
        } 
      }
    }  
     
    forall_in_edges(e,v)
    { if (flow[e] > 0)               // residual cap(e) > 0 
      { node w = G.source(e);
        PT  dp = pot[w];
      
        if (-cost[e] + pot[v] < dp)  // reduced cost cp(w,v) < 0 
        { cur[v] = e;
          STOP_TIMER(seek_time);
          return false;
        }
   
        dp -= -cost[e];  
 
        if (dp > cur_p) 
        { cur_p = dp;
          cur_e = e;
        } 
      }
    }  
    
    if (cur_p > min_p)
    { cur[v] = cur_e;
      pot[v] = cur_p - epsilon;      
    }
    else
    { if (excess[v] == 0)
        pot[v] = min_p; 
      else
      { cerr << "mincostflow::relabel(): problem is not feasible!\n";
        exit(1);
      }
    }  
    
    num_relabels[num_refines]++;       
    STOP_TIMER(relabel_time);

    return true;
  }

  
  // ---------------------------------------------------------
  // discharge:
  // preforms push and relabel operation on v as long as 
  // v is inactive, i. e. excess[v] == 0
  // ---------------------------------------------------------
  
  template <class cost_array, class flow_array, class node_slist>
  void discharge(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, node v)
  {    
    START_TIMER;
    
    if (!is_admissible(G,cost,cap,flow,pot,cur,v))
      relabel(G,cost,cap,flow,excess,pot,cur,v);
    
    while (excess[v] > 0)
    { edge e = cur[v];
      node w = G.opposite(e,v);
   
      if (excess[w] >= 0)
      { if (is_admissible(G,cost,cap,flow,pot,cur,w) || !relabel(G,cost,cap,flow,excess,pot,cur,w))
          push(G,cap,flow,excess,active,v,w,e); 
        else 
          if (excess[w] > 0)
            push_back(G,cap,flow,excess,v,w,e);
      }
      else 
      { push(G,cap,flow,excess,active,v,w,e);        
        if (excess[w] > 0 && !is_admissible(G,cost,cap,flow,pot,cur,w)) 
          relabel(G,cost,cap,flow,excess,pot,cur,w);        
      }
      
      if (excess[v] == 0) break;
      
      relabel(G,cost,cap,flow,excess,pot,cur,v);    
    }   
   
    num_discharges[num_refines]++;      
    STOP_TIMER(discharge_time);
  }  
  
  
  void reset_counter()
  {
    num_refines = 0;
    
    for(int i = 0; i < 16; i++) 
    { num_pushes[i] = 0;
      num_relabels[i] = 0;
      num_discharges[i] = 0;
    }
   
    total_flow = 0;
    total_cost = 0;
  
    num_nodes = 0;
    num_edges = 0;    

    init_time      = 0;
    refine_time    = 0;    
    seek_time      = 0;
    check_time     = 0;
    push_time      = 0;
    relabel_time   = 0;
    discharge_time = 0;
    total_time     = 0;      
  }
  
  
  template <class cost_array, class lcap_array, class ucap_array, class sup_array>
  void init(const graph& G, const cost_array& cost, const lcap_array& lcap, 
                            const ucap_array& ucap, const sup_array& sup, 
                                  cap_array&  cap, exc_array& exc, cur_array& cur)
  { START_TIMER;
      
    reset_counter();
      
    num_nodes = G.number_of_nodes();
    num_edges = G.number_of_edges();

    exc.use_node_data(G);  
    cur.use_node_data(G);

    node v;
    forall_nodes(v,G)
    { exc[v] = sup[v];
      cur[v] = G.outdeg(v) > 0 ? (edge) First_Adj_Edge(v,0) : (edge) First_Adj_Edge(v,1);
    }

    NT C = 0;

    cap.use_edge_data(G); 
  
    edge e;
    forall_edges(e,G)
    { cap[e] = ucap[e];  
  
      NT bound = lcap[e];
      if (bound != 0)
      { node v = G.source(e); 
        node w = G.target(e); 
        exc[v] -= bound;
        exc[w] += bound;
        cap[e] -= bound;            
      }
       
      if (cost[e] > C) C = cost[e];           
    } 

    epsilon = 1;
    while (epsilon < C) epsilon *= 2; 
    
    STOP_TIMER(init_time);
  }
};  

LEDA_END_NAMESPACE

#endif