feasible_flow1.t

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

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

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

#if !defined(LEDA_ROOT_INCL_ID)
#define LEDA_ROOT_INCL_ID 450447
#include <LEDA/PREAMBLE.h>
#endif

//-----------------------------------------------------------------------------
//
// FEASIBLE FLOW
//
// preflow-push + highest level
//
//-----------------------------------------------------------------------------

#include <LEDA/std/assert.h>


LEDA_BEGIN_NAMESPACE

#define NULL_NODE  ((node)0xFFFFFFFF)

template <class graph, class succ_array>
class ff1_node_level_queue {

typedef typename graph::node node;
typedef typename graph::edge edge;

node*  T;
node*  head;
node*  max;

succ_array& NEXT;

public:



ff1_node_level_queue(const graph& G, succ_array& succ) : NEXT(succ)
{ 
  int sz = G.number_of_nodes();
  T = new node[sz+1];
  for(node* p = T+sz; p >= T;  p--) *p = nil; 
  head = T+1;
  max = T;
  T[0] = NULL_NODE; 
 }

~ff1_node_level_queue() { delete[] T; }

void set_max(int i) { max = head + i; }

void insert(node v, int i)
{ // insert on level i
  node* h = head + i;
  if (h > max) max = h;
  NEXT[v] = *h;
  *h = v;
 }


void insert_non_max(node v, int i)
{ // do not update max pointer
  node* h = head + i;
  NEXT[v] = *h;
  *h = v;
 }


void insert_max(node v, int i)
{ // new max at level i
  max = head + i;
  *max = v;
  NEXT[v] = 0;
 }

node find_max()
{ while (*max == nil) max--;
  return *max;
}

node del_max()
{ node v = *max;
  *max = NEXT[v];
  return v; 
}

void del_max(node v) { *max = NEXT[v]; }


void clear() { while (max >= head) *max-- =  nil; }


};



template<class NT, class graph,
#if !defined(_MSC_VER)
         class succ_array   = node_array<typename graph::node,graph>,
#else
         class succ_array   = node_array<graph::node,graph>,
#endif
         class excess_array = node_array<int,graph>,
         class dist_array   = node_array<int,graph> >
class feasible_flow1  
{
  const int node_slots;
  const int edge_slots;

  typedef typename graph::node node;
  typedef typename graph::edge edge;


  float h;

  // statistics

  int num_pushes;
  int num_relabels;
  int num_updates;
  int num_inspections;

  float cputime;

  node source(const graph& G, edge e, node t)
  { if (graph::category() == opposite_graph_category)
      return G.opposite(e,t);
    else
      return G.source(e);
  }
  
  node target(const graph& G, edge e, node s)
  { if (graph::category() == opposite_graph_category)
      return G.opposite(e,s);
    else
      return G.target(e);
  }

public:

feasible_flow1() : node_slots(3), edge_slots(0), h(5.0)
{ reset_counters(); }

int pushes()      { return num_pushes; }
int relabels()    { return num_relabels; }
int updates()     { return num_updates; }
int inspections() { return num_inspections; }
float cpu_time()  { return cputime; }

void set_heuristic(float hh) { h = hh; }


void reset_counters()
{ num_pushes = 0;
  num_relabels = 0;
  num_updates = 0;
  num_inspections = 0;
  cputime = 0;
}




void compute_dist0(const graph& G, 
                   node t,
                   excess_array& excess,
                   succ_array& succ,
                   dist_array& dist,
                   ff1_node_level_queue<graph,succ_array>& U)
{ 
  // compute exact distance values by a "backward" bfs in the
  // residual network starting at t

  num_updates++;

  int n = G.number_of_nodes();

  dist.init(G,n);


  node level_last  = 0;
  node queue_last  = 0;

  dist[t] = 0;
  queue_last = t;
  level_last = t;

  int d = 0;

  node v = t;

  for(;;)
  { 
    edge e;
    forall_in_edges(e,v)
    { node u = source(G,e,v); 
      int& du = dist[u];
      if (du == n)
      { du = d+1;
        succ[queue_last] = u;
        queue_last = u;
       }
    }

    node next = succ[v];

    if (excess[v] > 0) U.insert(v,d);

    if (v == level_last)             
    { // finish level d
      if (queue_last == level_last) break; // no level d+1: stop
      level_last = queue_last;
     }
    v = next;
  }
}




template<class cap_array, class flow_array, class n_level_queue>
void compute_dist1(const graph& G, 
                   node t,
                   const cap_array& cap,
                   const flow_array& flow,
                   excess_array& excess,
                   succ_array& succ,
                   dist_array& dist,
                   n_level_queue& U)
{ 
  // compute exact distance values by a "backward" bfs in the
  // residual network starting at demand nodes

  num_updates++;

  int n = G.number_of_nodes();
  dist.init(G,n);

  node level_last  = 0;
  node queue_last  = 0;

  dist[t] = 0;
  queue_last = t;
  level_last = t;

  int d = 0;

  node v = t;

  for(;;)
  { 
    edge e;
    forall_out_edges(e,v)
    { if (flow[e] == 0) continue;
      node u = target(G,e,v); 
      int& du = dist[u];
      if (du == n)
      { du = d+1;
        succ[queue_last] = u;
        queue_last = u;
       }
     }

    forall_in_edges(e,v)
    { if (cap[e] == flow[e]) continue;
      node u = source(G,e,v); 
      int& du = dist[u];
      if (du == n)
      { du = d+1;
        succ[queue_last] = u;
        queue_last = u;
       }
    }

    node next = succ[v];

    if (excess[v] > 0) U.insert(v,d);

    if (v == level_last)             
    { // finish level d
      if (queue_last == level_last) break; // no level d+1: stop
      level_last = queue_last;
     }
    v = next;
  }
}



template<class cap_array, class flow_array>
bool run(const graph& G, node s, node t, int F,
                         const cap_array& cap, 
                         flow_array& flow)
{ 
  
  float T = used_time();

  reset_counters();

  flow.init(G,0);

  excess_array excess;
  excess.use_node_data(G,0);

  excess[s] =  F;
  excess[t] = -F;

  succ_array succ;
  succ.use_node_data(G);

  dist_array dist;
  dist.use_node_data(G);

  
  int n = G.number_of_nodes();
  int m = G.number_of_edges();

  int heuristic  = (h > 0) ? int(h*m) : int(-h*n);
  int limit_heur = heuristic;

  int relabel_count = 0;
  int push_count = 0;
  int insp_count = 0;

  ff1_node_level_queue<graph,succ_array>  U(G,succ);

  compute_dist0(G,t,excess,succ,dist,U);

  node v;

  while ((v = U.find_max()) != NULL_NODE) 
  {
    U.del_max(v);

    NT  ev = excess[v]; 
    int dv = dist[v]; 

    int  dmin = n;    
    edge emin = 0;
    NT   rmin = 0;

    edge e;
    forall_out_edges(e,v)
    { insp_count++;
      NT fe = flow[e];
      NT rc = cap[e] - fe;
      if (rc == 0) continue;
      node w = target(G,e,v);
      int dw = dist[w];
      if (dw == dv-1)
      { push_count++;
        NT ew = excess[w];
        if (ew == 0) U.insert_non_max(w,dw);
        if (ev <= rc) 
        { flow[e] = fe + ev; 
          excess[w] = ew + ev;
          excess[v] = 0; 
          goto NEXT_NODE;
         }
        flow[e] = fe + rc;
        excess[w] = ew + rc;
        ev -= rc;
        continue;
       }
      if (v != w && dw < dmin) { dmin = dw; emin = e; rmin = rc; }
    }

    forall_in_edges(e,v)
    { insp_count++;
      NT fe = flow[e];
      if (fe == 0) continue;
      node w = source(G,e,v);
      int dw = dist[w];
      if (dw == dv-1)
      { push_count++;
        NT ew = excess[w];
        if (ew == 0) U.insert_non_max(w,dw);
        if (ev <= fe) 
        { flow[e] = fe - ev; 
          excess[w] = ew + ev;
          excess[v] = 0; 
          goto NEXT_NODE;
         }
        flow[e] = 0;
        excess[w] = ew + fe;
        ev -= fe;
        continue;
       }
      if (v != w && dw < dmin) { dmin = dw; emin = e; rmin = -fe; }
    } 

    excess[v] = ev;

    // remaining excess at v
    // relabel vertex v (i.e. update dist[v]) because all
    // admissible edges in the residual graph have been saturated 

    relabel_count++;

    dv = dmin + 1;
    dist[v] = dv;

    if (dv >= n) goto NEXT_NODE;

    if (ev <= rmin || ev <= -rmin)
    { node w = G.opposite(emin,v);
      if (excess[w] == 0) U.insert(w,dmin); // possible new max
      excess[w] += ev;
      flow[emin] += (rmin > 0) ? ev : -ev;
      excess[v]  = 0;
      goto NEXT_NODE;
     }

    if ((h > 0 && insp_count    >= limit_heur) || 
        (h < 0 && relabel_count >= limit_heur))
    { 
      U.clear();
      compute_dist1(G,t,cap,flow,excess,succ,dist,U);
      limit_heur += heuristic;
      goto NEXT_NODE;
     }


   U.insert_max(v,dv);

NEXT_NODE:;
  }

  num_pushes = push_count;
  num_relabels = relabel_count;
  num_inspections = insp_count;

  cputime = used_time(T);

  return excess[t] == 0;
}





void statistics(ostream& out)
{ 
  out << string("%8d pushes %7d relabels %3d updates %9d insp  %.2f s",
                  pushes(), relabels(), updates(), inspections(), cpu_time());
  out << endl;

}


template<class excess, class cap_array, class flow_array>
bool check(const graph& G, const excess& excess,
                           const cap_array& cap, 
                           const flow_array& flow,
                           string& msg)
{
  edge e;
  forall_edges(e,G) 
  { if (flow[e] < 0 || flow[e] > cap[e]) 
    { msg = string("illegal flow value for edge %d",G.index(e));
      return false;
     }
   }

  node v;
  forall_nodes(v,G) 
  { edge e;
    forall_out_edges(e,v) 
    { node w = target(G,e,v);
      excess[v] -= flow[e]; 
      excess[w] += flow[e];
     }
  }

  int pos = 0;
  int neg = 0;

  forall_nodes(v,G)
  { if (excess[v] > 0) pos++;
    if (excess[v] < 0) neg++;
   }

  if (pos != 0 || neg != 0)
  { msg = string("non-zero excess: %d positive  %d negative",pos,neg);
    return false;
   }

  return true;
}


};




LEDA_END_NAMESPACE

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