src/work/jacint/max_flow.h
changeset 629 6620dfc606af
parent 602 580b329c2a0c
child 631 26819ef1611f
equal deleted inserted replaced
9:18f668399e53 10:47beb09fbab1
     1 // -*- C++ -*-
     1 // -*- C++ -*-
     2 
     2 
     3 /*
     3 /*
     4   Heuristics: 
     4   Heuristics:
     5   2 phase
     5   2 phase
     6   gap
     6   gap
     7   list 'level_list' on the nodes on level i implemented by hand
     7   list 'level_list' on the nodes on level i implemented by hand
     8   stack 'active' on the active nodes on level i
     8   stack 'active' on the active nodes on level i
     9   runs heuristic 'highest label' for H1*n relabels
     9   runs heuristic 'highest label' for H1*n relabels
    10   runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
    10   runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
    11  
    11 
    12   Parameters H0 and H1 are initialized to 20 and 1.
    12   Parameters H0 and H1 are initialized to 20 and 1.
    13 
    13 
    14   Constructors:
    14   Constructors:
    15 
    15 
    16   Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if 
    16   Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if
    17   FlowMap is not constant zero, and should be true if it is
    17   FlowMap is not constant zero, and should be true if it is
    18 
    18 
    19   Members:
    19   Members:
    20 
    20 
    21   void run()
    21   void run()
    22 
    22 
    23   Num flowValue() : returns the value of a maximum flow
    23   Num flowValue() : returns the value of a maximum flow
    24 
    24 
    25   void minMinCut(CutMap& M) : sets M to the characteristic vector of the 
    25   void minMinCut(CutMap& M) : sets M to the characteristic vector of the
    26   minimum min cut. M should be a map of bools initialized to false. ??Is it OK?
    26   minimum min cut. M should be a map of bools initialized to false. ??Is it OK?
    27 
    27 
    28   void maxMinCut(CutMap& M) : sets M to the characteristic vector of the 
    28   void maxMinCut(CutMap& M) : sets M to the characteristic vector of the
    29   maximum min cut. M should be a map of bools initialized to false.
    29   maximum min cut. M should be a map of bools initialized to false.
    30 
    30 
    31   void minCut(CutMap& M) : sets M to the characteristic vector of 
    31   void minCut(CutMap& M) : sets M to the characteristic vector of
    32   a min cut. M should be a map of bools initialized to false.
    32   a min cut. M should be a map of bools initialized to false.
    33 
    33 
    34 */
    34 */
    35 
    35 
    36 #ifndef HUGO_MAX_FLOW_H
    36 #ifndef HUGO_MAX_FLOW_H
    37 #define HUGO_MAX_FLOW_H
    37 #define HUGO_MAX_FLOW_H
    38 
       
    39 #define H0 20
       
    40 #define H1 1
       
    41 
    38 
    42 #include <vector>
    39 #include <vector>
    43 #include <queue>
    40 #include <queue>
    44 #include <stack>
    41 #include <stack>
    45 
    42 
    48 #include <hugo/invalid.h>
    45 #include <hugo/invalid.h>
    49 #include <hugo/maps.h>
    46 #include <hugo/maps.h>
    50 #include <for_each_macros.h>
    47 #include <for_each_macros.h>
    51 
    48 
    52 /// \file
    49 /// \file
    53 /// \brief Dimacs file format reader.
    50 /// \brief Maximum flows.
       
    51 /// \ingroup galgs
    54 
    52 
    55 namespace hugo {
    53 namespace hugo {
    56 
    54 
    57 
    55 
    58   //  ///\author Marton Makai, Jacint Szabo
    56   //  ///\author Marton Makai, Jacint Szabo
    59   /// A class for computing max flows and related quantities.
    57   /// A class for computing max flows and related quantities.
    60   template <typename Graph, typename Num, 
    58   /// \ingroup galgs
    61 	    typename CapMap=typename Graph::template EdgeMap<Num>, 
    59   template <typename Graph, typename Num,
       
    60 	    typename CapMap=typename Graph::template EdgeMap<Num>,
    62             typename FlowMap=typename Graph::template EdgeMap<Num> >
    61             typename FlowMap=typename Graph::template EdgeMap<Num> >
    63   class MaxFlow {
    62   class MaxFlow {
    64     
    63   protected:
    65     typedef typename Graph::Node Node;
    64     typedef typename Graph::Node Node;
    66     typedef typename Graph::NodeIt NodeIt;
    65     typedef typename Graph::NodeIt NodeIt;
    67     typedef typename Graph::OutEdgeIt OutEdgeIt;
    66     typedef typename Graph::OutEdgeIt OutEdgeIt;
    68     typedef typename Graph::InEdgeIt InEdgeIt;
    67     typedef typename Graph::InEdgeIt InEdgeIt;
    69 
    68 
    72     typedef typename std::vector<Node> VecNode;
    71     typedef typename std::vector<Node> VecNode;
    73 
    72 
    74     const Graph* g;
    73     const Graph* g;
    75     Node s;
    74     Node s;
    76     Node t;
    75     Node t;
    77     const CapMap* capacity;  
    76     const CapMap* capacity;
    78     FlowMap* flow;
    77     FlowMap* flow;
    79     int n;      //the number of nodes of G
    78     int n;      //the number of nodes of G
    80     typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
    79     typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
    81     typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
    80     typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
    82     typedef typename ResGW::Edge ResGWEdge;
    81     typedef typename ResGW::Edge ResGWEdge;
    83     //typedef typename ResGW::template NodeMap<bool> ReachedMap;
    82     //typedef typename ResGW::template NodeMap<bool> ReachedMap;
    84     typedef typename Graph::template NodeMap<int> ReachedMap;
    83     typedef typename Graph::template NodeMap<int> ReachedMap;
    85     ReachedMap level;
    84     ReachedMap level;
    86     //level works as a bool map in augmenting path algorithms 
    85     //level works as a bool map in augmenting path algorithms
    87     //and is used by bfs for storing reached information.
    86     //and is used by bfs for storing reached information.
    88     //In preflow, it shows levels of nodes.
    87     //In preflow, it shows levels of nodes.
    89     //typename Graph::template NodeMap<int> level;    
    88     //typename Graph::template NodeMap<int> level;
    90     typename Graph::template NodeMap<Num> excess; 
    89     typename Graph::template NodeMap<Num> excess;
    91     //   protected:
    90     //   protected:
    92     //     MaxFlow() { }
    91     //     MaxFlow() { }
    93     //     void set(const Graph& _G, Node _s, Node _t, const CapMap& _capacity, 
    92     //     void set(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,
    94     // 	     FlowMap& _flow) 
    93     // 	     FlowMap& _flow)
    95     //       {
    94     //       {
    96     // 	g=&_G; 
    95     // 	g=&_G;
    97     // 	s=_s; 
    96     // 	s=_s;
    98     // 	t=_t; 
    97     // 	t=_t;
    99     // 	capacity=&_capacity;
    98     // 	capacity=&_capacity;
   100     // 	flow=&_flow;
    99     // 	flow=&_flow;
   101     // 	n=_G.nodeNum;
   100     // 	n=_G.nodeNum;
   102     // 	level.set (_G); //kellene vmi ilyesmi fv 
   101     // 	level.set (_G); //kellene vmi ilyesmi fv
   103     // 	excess(_G,0); //itt is
   102     // 	excess(_G,0); //itt is
   104     //       }
   103     //       }
   105 
   104 
       
   105     // constants used for heuristics
       
   106     static const int H0=20;
       
   107     static const int H1=1;
       
   108 
   106   public:
   109   public:
   107  
   110 
   108     ///\todo Document this.
   111     ///\todo Document this.
   109     ///\todo Maybe, it should be PRE_FLOW instead.
   112     ///\todo Maybe, it should be PRE_FLOW instead.
       
   113     ///- \c NO_FLOW means nothing,
   110     ///- \c ZERO_FLOW means something,
   114     ///- \c ZERO_FLOW means something,
   111     ///- \c GEN_FLOW means something else,
   115     ///- \c GEN_FLOW means something else,
   112     ///- \c PREFLOW is something different.
   116     ///- \c PRE_FLOW is something different.
   113     enum flowEnum{
   117     enum flowEnum{
   114       ZERO_FLOW=0,
   118       ZERO_FLOW,
   115       GEN_FLOW=1,
   119       GEN_FLOW,
   116       PREFLOW=2
   120       PRE_FLOW,
       
   121       NO_FLOW
   117     };
   122     };
   118 
   123 
   119     MaxFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity, 
   124     MaxFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,
   120 	    FlowMap& _flow) :
   125 	    FlowMap& _flow) :
   121       g(&_G), s(_s), t(_t), capacity(&_capacity), 
   126       g(&_G), s(_s), t(_t), capacity(&_capacity),
   122       flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {}
   127       flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {}
   123 
   128 
   124     /// A max flow algorithm is run.
   129     /// A max flow algorithm is run.
   125     ///\pre the flow have to be 0 at the beginning.
   130     /// \pre The flow have to satisfy the requirements
   126     void run() {
   131     /// stated in fe.
   127       preflow(ZERO_FLOW);
   132     void run(flowEnum fe=ZERO_FLOW) {
   128     }
   133       preflow(fe);
   129     
   134     }
   130     /// A preflow algorithm is run. 
   135 
   131     ///\pre The initial edge-map have to be a 
   136     /// A preflow algorithm is run.
       
   137     /// \pre The initial edge-map have to be a
   132     /// zero flow if \c fe is \c ZERO_FLOW,
   138     /// zero flow if \c fe is \c ZERO_FLOW,
   133     /// a flow if \c fe is \c GEN_FLOW, 
   139     /// a flow if \c fe is \c GEN_FLOW,
   134     /// and a pre-flow it is \c PREFLOW.
   140     /// a pre-flow if fe is \c PRE_FLOW and
       
   141     /// anything if fe is NO_FLOW.
   135     void preflow(flowEnum fe) {
   142     void preflow(flowEnum fe) {
   136       preflowPhase0(fe);
   143       preflowPhase0(fe);
   137       preflowPhase1();
   144       preflowPhase1();
   138     }
   145     }
   139 
   146 
   140     /// Run the first phase of preflow, starting from a 0 flow, from a flow, 
   147     /// Run the first phase of preflow, starting from a 0 flow, from a flow,
   141     /// or from a preflow, according to \c fe.
   148     /// or from a preflow, of from undefined value according to \c fe.
   142     void preflowPhase0( flowEnum fe );
   149     void preflowPhase0(flowEnum fe);
   143 
   150 
   144     /// Second phase of preflow.
   151     /// Second phase of preflow.
   145     void preflowPhase1();
   152     void preflowPhase1();
   146 
   153 
   147     /// Starting from a flow, this method searches for an augmenting path 
   154     /// Starting from a flow, this method searches for an augmenting path
   148     /// according to the Edmonds-Karp algorithm 
   155     /// according to the Edmonds-Karp algorithm
   149     /// and augments the flow on if any. 
   156     /// and augments the flow on if any.
   150     /// The return value shows if the augmentation was succesful.
   157     /// The return value shows if the augmentation was succesful.
   151     bool augmentOnShortestPath();
   158     bool augmentOnShortestPath();
   152 
   159 
   153     /// Starting from a flow, this method searches for an augmenting blockin 
   160     /// Starting from a flow, this method searches for an augmenting blocking
   154     /// flow according to Dinits' algorithm and augments the flow on if any. 
   161     /// flow according to Dinits' algorithm and augments the flow on if any.
   155     /// The blocking flow is computed in a physically constructed 
   162     /// The blocking flow is computed in a physically constructed
   156     /// residual graph of type \c Mutablegraph.
   163     /// residual graph of type \c Mutablegraph.
   157     /// The return value show sif the augmentation was succesful.
   164     /// The return value show sif the augmentation was succesful.
   158     template<typename MutableGraph> bool augmentOnBlockingFlow();
   165     template<typename MutableGraph> bool augmentOnBlockingFlow();
   159 
   166 
   160     /// The same as \c augmentOnBlockingFlow<MutableGraph> but the 
   167     /// The same as \c augmentOnBlockingFlow<MutableGraph> but the
   161     /// residual graph is not constructed physically.
   168     /// residual graph is not constructed physically.
   162     /// The return value shows if the augmentation was succesful.
   169     /// The return value shows if the augmentation was succesful.
   163     bool augmentOnBlockingFlow2();
   170     bool augmentOnBlockingFlow2();
   164 
   171 
   165     /// Returns the actual flow value.
   172     /// Returns the actual flow value.
   166     /// More precisely, it returns the negative excess of s, thus 
   173     /// More precisely, it returns the negative excess of s, thus
   167     /// this works also for preflows.
   174     /// this works also for preflows.
   168     Num flowValue() { 
   175     Num flowValue() {
   169       Num a=0;
   176       Num a=0;
   170       FOR_EACH_INC_LOC(OutEdgeIt, e, *g, s) a+=(*flow)[e];
   177       FOR_EACH_INC_LOC(InEdgeIt, e, *g, t) a+=(*flow)[e];
   171       FOR_EACH_INC_LOC(InEdgeIt, e, *g, s) a-=(*flow)[e];
   178       FOR_EACH_INC_LOC(OutEdgeIt, e, *g, t) a-=(*flow)[e];
   172       return a;
   179       return a;
   173     }
   180     }
   174 
   181 
   175     /// Should be used between preflowPhase0 and preflowPhase1.
   182     /// Should be used between preflowPhase0 and preflowPhase1.
   176     ///\todo We have to make some status variable which shows the actual state 
   183     /// \todo We have to make some status variable which shows the
   177     /// of the class. This enables us to determine which methods are valid 
   184     /// actual state
       
   185     /// of the class. This enables us to determine which methods are valid
   178     /// for MinCut computation
   186     /// for MinCut computation
   179     template<typename _CutMap>
   187     template<typename _CutMap>
   180     void actMinCut(_CutMap& M) {
   188     void actMinCut(_CutMap& M) {
   181       NodeIt v;
   189       NodeIt v;
   182       for(g->first(v); g->valid(v); g->next(v)) {
   190       for(g->first(v); g->valid(v); g->next(v)) {
   186 	  M.set(v,true);
   194 	  M.set(v,true);
   187 	}
   195 	}
   188       }
   196       }
   189     }
   197     }
   190 
   198 
   191     /// The unique inclusionwise minimum cut is computed by 
   199     /// The unique inclusionwise minimum cut is computed by
   192     /// processing a bfs from s in the residual graph.
   200     /// processing a bfs from s in the residual graph.
   193     ///\pre flow have to be a max flow otherwise it will the whole node-set.
   201     /// \pre flow have to be a max flow otherwise it will the whole node-set.
   194     template<typename _CutMap>
   202     template<typename _CutMap>
   195     void minMinCut(_CutMap& M) {
   203     void minMinCut(_CutMap& M) {
   196     
   204 
   197       std::queue<Node> queue;
   205       std::queue<Node> queue;
   198       
   206 
   199       M.set(s,true);      
   207       M.set(s,true);
   200       queue.push(s);
   208       queue.push(s);
   201 
   209 
   202       while (!queue.empty()) {
   210       while (!queue.empty()) {
   203         Node w=queue.front();
   211         Node w=queue.front();
   204 	queue.pop();
   212 	queue.pop();
   208 	  Node v=g->head(e);
   216 	  Node v=g->head(e);
   209 	  if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
   217 	  if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
   210 	    queue.push(v);
   218 	    queue.push(v);
   211 	    M.set(v, true);
   219 	    M.set(v, true);
   212 	  }
   220 	  }
   213 	} 
   221 	}
   214 
   222 
   215 	InEdgeIt f;
   223 	InEdgeIt f;
   216 	for(g->first(f,w) ; g->valid(f); g->next(f)) {
   224 	for(g->first(f,w) ; g->valid(f); g->next(f)) {
   217 	  Node v=g->tail(f);
   225 	  Node v=g->tail(f);
   218 	  if (!M[v] && (*flow)[f] > 0 ) {
   226 	  if (!M[v] && (*flow)[f] > 0 ) {
   219 	    queue.push(v);
   227 	    queue.push(v);
   220 	    M.set(v, true);
   228 	    M.set(v, true);
   221 	  }
   229 	  }
   222 	} 
   230 	}
   223       }
   231       }
   224     }
   232     }
   225 
   233 
   226 
   234 
   227     /// The unique inclusionwise maximum cut is computed by 
   235     /// The unique inclusionwise maximum cut is computed by
   228     /// processing a reverse bfs from t in the residual graph.
   236     /// processing a reverse bfs from t in the residual graph.
   229     ///\pre flow have to be a max flow otherwise it will be empty.
   237     /// \pre flow have to be a max flow otherwise it will be empty.
   230     template<typename _CutMap>
   238     template<typename _CutMap>
   231     void maxMinCut(_CutMap& M) {
   239     void maxMinCut(_CutMap& M) {
   232 
   240 
   233       NodeIt v;
   241       NodeIt v;
   234       for(g->first(v) ; g->valid(v); g->next(v)) {
   242       for(g->first(v) ; g->valid(v); g->next(v)) {
   235 	M.set(v, true);
   243 	M.set(v, true);
   236       }
   244       }
   237 
   245 
   238       std::queue<Node> queue;
   246       std::queue<Node> queue;
   239       
   247 
   240       M.set(t,false);        
   248       M.set(t,false);
   241       queue.push(t);
   249       queue.push(t);
   242 
   250 
   243       while (!queue.empty()) {
   251       while (!queue.empty()) {
   244         Node w=queue.front();
   252         Node w=queue.front();
   245 	queue.pop();
   253 	queue.pop();
   246 
       
   247 
   254 
   248 	InEdgeIt e;
   255 	InEdgeIt e;
   249 	for(g->first(e,w) ; g->valid(e); g->next(e)) {
   256 	for(g->first(e,w) ; g->valid(e); g->next(e)) {
   250 	  Node v=g->tail(e);
   257 	  Node v=g->tail(e);
   251 	  if (M[v] && (*flow)[e] < (*capacity)[e] ) {
   258 	  if (M[v] && (*flow)[e] < (*capacity)[e] ) {
   252 	    queue.push(v);
   259 	    queue.push(v);
   253 	    M.set(v, false);
   260 	    M.set(v, false);
   254 	  }
   261 	  }
   255 	}
   262 	}
   256 	
   263 
   257 	OutEdgeIt f;
   264 	OutEdgeIt f;
   258 	for(g->first(f,w) ; g->valid(f); g->next(f)) {
   265 	for(g->first(f,w) ; g->valid(f); g->next(f)) {
   259 	  Node v=g->head(f);
   266 	  Node v=g->head(f);
   260 	  if (M[v] && (*flow)[f] > 0 ) {
   267 	  if (M[v] && (*flow)[f] > 0 ) {
   261 	    queue.push(v);
   268 	    queue.push(v);
   272 
   279 
   273     ///
   280     ///
   274     void resetSource(Node _s) { s=_s; }
   281     void resetSource(Node _s) { s=_s; }
   275     ///
   282     ///
   276     void resetTarget(Node _t) { t=_t; }
   283     void resetTarget(Node _t) { t=_t; }
   277    
   284 
   278     /// capacity-map is changed.
   285     /// capacity-map is changed.
   279     void resetCap(const CapMap& _cap) { capacity=&_cap; }
   286     void resetCap(const CapMap& _cap) { capacity=&_cap; }
   280     
   287 
   281     /// flow-map is changed. 
   288     /// flow-map is changed.
   282     void resetFlow(FlowMap& _flow) { flow=&_flow; }
   289     void resetFlow(FlowMap& _flow) { flow=&_flow; }
   283 
   290 
   284 
   291 
   285   private:
   292   private:
   286 
   293 
   287     int push(Node w, VecStack& active) {
   294     int push(Node w, VecStack& active) {
   288       
   295 
   289       int lev=level[w];
   296       int lev=level[w];
   290       Num exc=excess[w];
   297       Num exc=excess[w];
   291       int newlevel=n;       //bound on the next level of w
   298       int newlevel=n;       //bound on the next level of w
   292 	  
   299 
   293       OutEdgeIt e;
   300       OutEdgeIt e;
   294       for(g->first(e,w); g->valid(e); g->next(e)) {
   301       for(g->first(e,w); g->valid(e); g->next(e)) {
   295 	    
   302 
   296 	if ( (*flow)[e] >= (*capacity)[e] ) continue; 
   303 	if ( (*flow)[e] >= (*capacity)[e] ) continue;
   297 	Node v=g->head(e);            
   304 	Node v=g->head(e);
   298 	    
   305 
   299 	if( lev > level[v] ) { //Push is allowed now
   306 	if( lev > level[v] ) { //Push is allowed now
   300 	  
   307 
   301 	  if ( excess[v]<=0 && v!=t && v!=s ) {
   308 	  if ( excess[v]<=0 && v!=t && v!=s ) {
   302 	    int lev_v=level[v];
   309 	    int lev_v=level[v];
   303 	    active[lev_v].push(v);
   310 	    active[lev_v].push(v);
   304 	  }
   311 	  }
   305 	  
   312 
   306 	  Num cap=(*capacity)[e];
   313 	  Num cap=(*capacity)[e];
   307 	  Num flo=(*flow)[e];
   314 	  Num flo=(*flow)[e];
   308 	  Num remcap=cap-flo;
   315 	  Num remcap=cap-flo;
   309 	  
   316 
   310 	  if ( remcap >= exc ) { //A nonsaturating push.
   317 	  if ( remcap >= exc ) { //A nonsaturating push.
   311 	    
   318 
   312 	    flow->set(e, flo+exc);
   319 	    flow->set(e, flo+exc);
   313 	    excess.set(v, excess[v]+exc);
   320 	    excess.set(v, excess[v]+exc);
   314 	    exc=0;
   321 	    exc=0;
   315 	    break; 
   322 	    break;
   316 	    
   323 
   317 	  } else { //A saturating push.
   324 	  } else { //A saturating push.
   318 	    flow->set(e, cap);
   325 	    flow->set(e, cap);
   319 	    excess.set(v, excess[v]+remcap);
   326 	    excess.set(v, excess[v]+remcap);
   320 	    exc-=remcap;
   327 	    exc-=remcap;
   321 	  }
   328 	  }
   322 	} else if ( newlevel > level[v] ) newlevel = level[v];
   329 	} else if ( newlevel > level[v] ) newlevel = level[v];
   323       } //for out edges wv 
   330       } //for out edges wv
   324       
   331 
   325       if ( exc > 0 ) {	
   332       if ( exc > 0 ) {
   326 	InEdgeIt e;
   333 	InEdgeIt e;
   327 	for(g->first(e,w); g->valid(e); g->next(e)) {
   334 	for(g->first(e,w); g->valid(e); g->next(e)) {
   328 	  
   335 
   329 	  if( (*flow)[e] <= 0 ) continue; 
   336 	  if( (*flow)[e] <= 0 ) continue;
   330 	  Node v=g->tail(e); 
   337 	  Node v=g->tail(e);
   331 	  
   338 
   332 	  if( lev > level[v] ) { //Push is allowed now
   339 	  if( lev > level[v] ) { //Push is allowed now
   333 	    
   340 
   334 	    if ( excess[v]<=0 && v!=t && v!=s ) {
   341 	    if ( excess[v]<=0 && v!=t && v!=s ) {
   335 	      int lev_v=level[v];
   342 	      int lev_v=level[v];
   336 	      active[lev_v].push(v);
   343 	      active[lev_v].push(v);
   337 	    }
   344 	    }
   338 	    
   345 
   339 	    Num flo=(*flow)[e];
   346 	    Num flo=(*flow)[e];
   340 	    
   347 
   341 	    if ( flo >= exc ) { //A nonsaturating push.
   348 	    if ( flo >= exc ) { //A nonsaturating push.
   342 	      
   349 
   343 	      flow->set(e, flo-exc);
   350 	      flow->set(e, flo-exc);
   344 	      excess.set(v, excess[v]+exc);
   351 	      excess.set(v, excess[v]+exc);
   345 	      exc=0;
   352 	      exc=0;
   346 	      break; 
   353 	      break;
   347 	    } else {  //A saturating push.
   354 	    } else {  //A saturating push.
   348 	      
   355 
   349 	      excess.set(v, excess[v]+flo);
   356 	      excess.set(v, excess[v]+flo);
   350 	      exc-=flo;
   357 	      exc-=flo;
   351 	      flow->set(e,0);
   358 	      flow->set(e,0);
   352 	    }  
   359 	    }
   353 	  } else if ( newlevel > level[v] ) newlevel = level[v];
   360 	  } else if ( newlevel > level[v] ) newlevel = level[v];
   354 	} //for in edges vw
   361 	} //for in edges vw
   355 	
   362 
   356       } // if w still has excess after the out edge for cycle
   363       } // if w still has excess after the out edge for cycle
   357       
   364 
   358       excess.set(w, exc);
   365       excess.set(w, exc);
   359       
   366 
   360       return newlevel;
   367       return newlevel;
   361     }
   368     }
   362 
   369 
   363 
   370 
   364     void preflowPreproc ( flowEnum fe, VecStack& active, 
   371     void preflowPreproc(flowEnum fe, VecStack& active,
   365 			  VecNode& level_list, NNMap& left, NNMap& right ) 
   372 			VecNode& level_list, NNMap& left, NNMap& right)
   366     {
   373     {
   367 
       
   368       std::queue<Node> bfs_queue;
   374       std::queue<Node> bfs_queue;
   369       
   375 
   370       switch ( fe ) {
   376       switch (fe) {
   371       case ZERO_FLOW: 
   377       case ZERO_FLOW:
   372 	{
   378 	{
   373 	  //Reverse_bfs from t, to find the starting level.
   379 	  //Reverse_bfs from t, to find the starting level.
   374 	  level.set(t,0);
   380 	  level.set(t,0);
   375 	  bfs_queue.push(t);
   381 	  bfs_queue.push(t);
   376 	
   382 
   377 	  while (!bfs_queue.empty()) {
   383 	  while (!bfs_queue.empty()) {
   378 	    
   384 
   379 	    Node v=bfs_queue.front();	
   385 	    Node v=bfs_queue.front();
   380 	    bfs_queue.pop();
   386 	    bfs_queue.pop();
   381 	    int l=level[v]+1;
   387 	    int l=level[v]+1;
   382 	    
   388 
   383 	    InEdgeIt e;
   389 	    InEdgeIt e;
   384 	    for(g->first(e,v); g->valid(e); g->next(e)) {
   390 	    for(g->first(e,v); g->valid(e); g->next(e)) {
   385 	      Node w=g->tail(e);
   391 	      Node w=g->tail(e);
   386 	      if ( level[w] == n && w != s ) {
   392 	      if ( level[w] == n && w != s ) {
   387 		bfs_queue.push(w);
   393 		bfs_queue.push(w);
   391 		level_list[l]=w;
   397 		level_list[l]=w;
   392 		level.set(w, l);
   398 		level.set(w, l);
   393 	      }
   399 	      }
   394 	    }
   400 	    }
   395 	  }
   401 	  }
   396 	  
   402 
   397 	  //the starting flow
   403 	  //the starting flow
   398 	  OutEdgeIt e;
   404 	  OutEdgeIt e;
   399 	  for(g->first(e,s); g->valid(e); g->next(e)) 
   405 	  for(g->first(e,s); g->valid(e); g->next(e))
   400 	    {
   406 	    {
   401 	      Num c=(*capacity)[e];
   407 	      Num c=(*capacity)[e];
   402 	      if ( c <= 0 ) continue;
   408 	      if ( c <= 0 ) continue;
   403 	      Node w=g->head(e);
   409 	      Node w=g->head(e);
   404 	      if ( level[w] < n ) {	  
   410 	      if ( level[w] < n ) {
   405 		if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
   411 		if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
   406 		flow->set(e, c); 
   412 		flow->set(e, c);
   407 		excess.set(w, excess[w]+c);
   413 		excess.set(w, excess[w]+c);
   408 	      }
   414 	      }
   409 	    }
   415 	    }
   410 	  break;
   416 	  break;
   411 	}
   417 	}
   412 	
   418 
   413       case GEN_FLOW:
   419       case GEN_FLOW:
   414       case PREFLOW:
   420       case PRE_FLOW:
   415 	{
   421 	{
   416 	  //Reverse_bfs from t in the residual graph, 
   422 	  //Reverse_bfs from t in the residual graph,
   417 	  //to find the starting level.
   423 	  //to find the starting level.
   418 	  level.set(t,0);
   424 	  level.set(t,0);
   419 	  bfs_queue.push(t);
   425 	  bfs_queue.push(t);
   420 	  
   426 
   421 	  while (!bfs_queue.empty()) {
   427 	  while (!bfs_queue.empty()) {
   422 	    
   428 
   423 	    Node v=bfs_queue.front();	
   429 	    Node v=bfs_queue.front();
   424 	    bfs_queue.pop();
   430 	    bfs_queue.pop();
   425 	    int l=level[v]+1;
   431 	    int l=level[v]+1;
   426 	    
   432 
   427 	    InEdgeIt e;
   433 	    InEdgeIt e;
   428 	    for(g->first(e,v); g->valid(e); g->next(e)) {
   434 	    for(g->first(e,v); g->valid(e); g->next(e)) {
   429 	      if ( (*capacity)[e] <= (*flow)[e] ) continue;
   435 	      if ( (*capacity)[e] <= (*flow)[e] ) continue;
   430 	      Node w=g->tail(e);
   436 	      Node w=g->tail(e);
   431 	      if ( level[w] == n && w != s ) {
   437 	      if ( level[w] == n && w != s ) {
   435 		right.set(w,first);
   441 		right.set(w,first);
   436 		level_list[l]=w;
   442 		level_list[l]=w;
   437 		level.set(w, l);
   443 		level.set(w, l);
   438 	      }
   444 	      }
   439 	    }
   445 	    }
   440 	    
   446 
   441 	    OutEdgeIt f;
   447 	    OutEdgeIt f;
   442 	    for(g->first(f,v); g->valid(f); g->next(f)) {
   448 	    for(g->first(f,v); g->valid(f); g->next(f)) {
   443 	      if ( 0 >= (*flow)[f] ) continue;
   449 	      if ( 0 >= (*flow)[f] ) continue;
   444 	      Node w=g->head(f);
   450 	      Node w=g->head(f);
   445 	      if ( level[w] == n && w != s ) {
   451 	      if ( level[w] == n && w != s ) {
   450 		level_list[l]=w;
   456 		level_list[l]=w;
   451 		level.set(w, l);
   457 		level.set(w, l);
   452 	      }
   458 	      }
   453 	    }
   459 	    }
   454 	  }
   460 	  }
   455 	  
   461 
   456 	  
   462 
   457 	  //the starting flow
   463 	  //the starting flow
   458 	  OutEdgeIt e;
   464 	  OutEdgeIt e;
   459 	  for(g->first(e,s); g->valid(e); g->next(e)) 
   465 	  for(g->first(e,s); g->valid(e); g->next(e))
   460 	    {
   466 	    {
   461 	      Num rem=(*capacity)[e]-(*flow)[e];
   467 	      Num rem=(*capacity)[e]-(*flow)[e];
   462 	      if ( rem <= 0 ) continue;
   468 	      if ( rem <= 0 ) continue;
   463 	      Node w=g->head(e);
   469 	      Node w=g->head(e);
   464 	      if ( level[w] < n ) {	  
   470 	      if ( level[w] < n ) {
   465 		if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
   471 		if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
   466 		flow->set(e, (*capacity)[e]); 
   472 		flow->set(e, (*capacity)[e]);
   467 		excess.set(w, excess[w]+rem);
   473 		excess.set(w, excess[w]+rem);
   468 	      }
   474 	      }
   469 	    }
   475 	    }
   470 	  
   476 
   471 	  InEdgeIt f;
   477 	  InEdgeIt f;
   472 	  for(g->first(f,s); g->valid(f); g->next(f)) 
   478 	  for(g->first(f,s); g->valid(f); g->next(f))
   473 	    {
   479 	    {
   474 	      if ( (*flow)[f] <= 0 ) continue;
   480 	      if ( (*flow)[f] <= 0 ) continue;
   475 	      Node w=g->tail(f);
   481 	      Node w=g->tail(f);
   476 	      if ( level[w] < n ) {	  
   482 	      if ( level[w] < n ) {
   477 		if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
   483 		if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
   478 		excess.set(w, excess[w]+(*flow)[f]);
   484 		excess.set(w, excess[w]+(*flow)[f]);
   479 		flow->set(f, 0); 
   485 		flow->set(f, 0);
   480 	      }
   486 	      }
   481 	    }  
   487 	    }
   482 	  break;
   488 	  break;
   483 	} //case PREFLOW
   489 	} //case PRE_FLOW
   484       }
   490       }
   485     } //preflowPreproc
   491     } //preflowPreproc
   486 
   492 
   487 
   493 
   488 
   494 
   489     void relabel(Node w, int newlevel, VecStack& active,  
   495     void relabel(Node w, int newlevel, VecStack& active,
   490 		 VecNode& level_list, NNMap& left, 
   496 		 VecNode& level_list, NNMap& left,
   491 		 NNMap& right, int& b, int& k, bool what_heur ) 
   497 		 NNMap& right, int& b, int& k, bool what_heur )
   492     {
   498     {
   493 
   499 
   494       Num lev=level[w];	
   500       Num lev=level[w];
   495       
   501 
   496       Node right_n=right[w];
   502       Node right_n=right[w];
   497       Node left_n=left[w];
   503       Node left_n=left[w];
   498       
   504 
   499       //unlacing starts
   505       //unlacing starts
   500       if ( g->valid(right_n) ) {
   506       if ( g->valid(right_n) ) {
   501 	if ( g->valid(left_n) ) {
   507 	if ( g->valid(left_n) ) {
   502 	  right.set(left_n, right_n);
   508 	  right.set(left_n, right_n);
   503 	  left.set(right_n, left_n);
   509 	  left.set(right_n, left_n);
   504 	} else {
   510 	} else {
   505 	  level_list[lev]=right_n;   
   511 	  level_list[lev]=right_n;
   506 	  left.set(right_n, INVALID);
   512 	  left.set(right_n, INVALID);
   507 	} 
   513 	}
   508       } else {
   514       } else {
   509 	if ( g->valid(left_n) ) {
   515 	if ( g->valid(left_n) ) {
   510 	  right.set(left_n, INVALID);
   516 	  right.set(left_n, INVALID);
   511 	} else { 
   517 	} else {
   512 	  level_list[lev]=INVALID;   
   518 	  level_list[lev]=INVALID;
   513 	} 
   519 	}
   514       } 
   520       }
   515       //unlacing ends
   521       //unlacing ends
   516 		
   522 
   517       if ( !g->valid(level_list[lev]) ) {
   523       if ( !g->valid(level_list[lev]) ) {
   518 	      
   524 
   519 	//gapping starts
   525 	//gapping starts
   520 	for (int i=lev; i!=k ; ) {
   526 	for (int i=lev; i!=k ; ) {
   521 	  Node v=level_list[++i];
   527 	  Node v=level_list[++i];
   522 	  while ( g->valid(v) ) {
   528 	  while ( g->valid(v) ) {
   523 	    level.set(v,n);
   529 	    level.set(v,n);
   526 	  level_list[i]=INVALID;
   532 	  level_list[i]=INVALID;
   527 	  if ( !what_heur ) {
   533 	  if ( !what_heur ) {
   528 	    while ( !active[i].empty() ) {
   534 	    while ( !active[i].empty() ) {
   529 	      active[i].pop();    //FIXME: ezt szebben kene
   535 	      active[i].pop();    //FIXME: ezt szebben kene
   530 	    }
   536 	    }
   531 	  }	     
   537 	  }
   532 	}
   538 	}
   533 	
   539 
   534 	level.set(w,n);
   540 	level.set(w,n);
   535 	b=lev-1;
   541 	b=lev-1;
   536 	k=b;
   542 	k=b;
   537 	//gapping ends
   543 	//gapping ends
   538 	
   544 
   539       } else {
   545       } else {
   540 	
   546 
   541 	if ( newlevel == n ) level.set(w,n); 
   547 	if ( newlevel == n ) level.set(w,n);
   542 	else {
   548 	else {
   543 	  level.set(w,++newlevel);
   549 	  level.set(w,++newlevel);
   544 	  active[newlevel].push(w);
   550 	  active[newlevel].push(w);
   545 	  if ( what_heur ) b=newlevel;
   551 	  if ( what_heur ) b=newlevel;
   546 	  if ( k < newlevel ) ++k;      //now k=newlevel
   552 	  if ( k < newlevel ) ++k;      //now k=newlevel
   549 	  right.set(w,first);
   555 	  right.set(w,first);
   550 	  left.set(w,INVALID);
   556 	  left.set(w,INVALID);
   551 	  level_list[newlevel]=w;
   557 	  level_list[newlevel]=w;
   552 	}
   558 	}
   553       }
   559       }
   554       
   560 
   555     } //relabel
   561     } //relabel
   556 
   562 
   557 
   563 
   558     template<typename MapGraphWrapper> 
   564     template<typename MapGraphWrapper>
   559     class DistanceMap {
   565     class DistanceMap {
   560     protected:
   566     protected:
   561       const MapGraphWrapper* g;
   567       const MapGraphWrapper* g;
   562       typename MapGraphWrapper::template NodeMap<int> dist; 
   568       typename MapGraphWrapper::template NodeMap<int> dist;
   563     public:
   569     public:
   564       DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { }
   570       DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { }
   565       void set(const typename MapGraphWrapper::Node& n, int a) { 
   571       void set(const typename MapGraphWrapper::Node& n, int a) {
   566 	dist.set(n, a); 
   572 	dist.set(n, a);
   567       }
   573       }
   568       int operator[](const typename MapGraphWrapper::Node& n) 
   574       int operator[](const typename MapGraphWrapper::Node& n)
   569       { return dist[n]; }
   575       { return dist[n]; }
   570       //       int get(const typename MapGraphWrapper::Node& n) const { 
   576       //       int get(const typename MapGraphWrapper::Node& n) const {
   571       // 	return dist[n]; }
   577       // 	return dist[n]; }
   572       //       bool get(const typename MapGraphWrapper::Edge& e) const { 
   578       //       bool get(const typename MapGraphWrapper::Edge& e) const {
   573       // 	return (dist.get(g->tail(e))<dist.get(g->head(e))); }
   579       // 	return (dist.get(g->tail(e))<dist.get(g->head(e))); }
   574       bool operator[](const typename MapGraphWrapper::Edge& e) const { 
   580       bool operator[](const typename MapGraphWrapper::Edge& e) const {
   575 	return (dist[g->tail(e)]<dist[g->head(e)]); 
   581 	return (dist[g->tail(e)]<dist[g->head(e)]);
   576       }
   582       }
   577     };
   583     };
   578     
   584 
   579   };
   585   };
   580 
   586 
   581 
   587 
   582   template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   588   template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   583   void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase0( flowEnum fe ) 
   589   void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase0( flowEnum fe )
   584   {
   590   {
   585       
   591 
   586     int heur0=(int)(H0*n);  //time while running 'bound decrease' 
   592     int heur0=(int)(H0*n);  //time while running 'bound decrease'
   587     int heur1=(int)(H1*n);  //time while running 'highest label'
   593     int heur1=(int)(H1*n);  //time while running 'highest label'
   588     int heur=heur1;         //starting time interval (#of relabels)
   594     int heur=heur1;         //starting time interval (#of relabels)
   589     int numrelabel=0;
   595     int numrelabel=0;
   590      
   596 
   591     bool what_heur=1;       
   597     bool what_heur=1;
   592     //It is 0 in case 'bound decrease' and 1 in case 'highest label'
   598     //It is 0 in case 'bound decrease' and 1 in case 'highest label'
   593 
   599 
   594     bool end=false;     
   600     bool end=false;
   595     //Needed for 'bound decrease', true means no active nodes are above bound b.
   601     //Needed for 'bound decrease', true means no active nodes are above bound
       
   602     //b.
   596 
   603 
   597     int k=n-2;  //bound on the highest level under n containing a node
   604     int k=n-2;  //bound on the highest level under n containing a node
   598     int b=k;    //bound on the highest level under n of an active node
   605     int b=k;    //bound on the highest level under n of an active node
   599       
   606 
   600     VecStack active(n);
   607     VecStack active(n);
   601       
   608 
   602     NNMap left(*g, INVALID);
   609     NNMap left(*g, INVALID);
   603     NNMap right(*g, INVALID);
   610     NNMap right(*g, INVALID);
   604     VecNode level_list(n,INVALID);
   611     VecNode level_list(n,INVALID);
   605     //List of the nodes in level i<n, set to n.
   612     //List of the nodes in level i<n, set to n.
   606 
   613 
   607     NodeIt v;
   614     NodeIt v;
   608     for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
   615     for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
   609     //setting each node to level n
   616     //setting each node to level n
   610       
   617 
   611     switch ( fe ) {
   618     switch (fe) {
   612     case PREFLOW:
   619     case PRE_FLOW:
   613       {
   620       {
   614 	//counting the excess
   621 	//computing the excess
   615 	NodeIt v;
   622 	NodeIt v;
   616 	for(g->first(v); g->valid(v); g->next(v)) {
   623 	for(g->first(v); g->valid(v); g->next(v)) {
   617 	  Num exc=0;
   624 	  Num exc=0;
   618 	  
   625 
   619 	  InEdgeIt e;
   626 	  InEdgeIt e;
   620 	  for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
   627 	  for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
   621 	  OutEdgeIt f;
   628 	  OutEdgeIt f;
   622 	  for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
   629 	  for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
   623 	    
   630 
   624 	  excess.set(v,exc);	  
   631 	  excess.set(v,exc);
   625 	    
   632 
   626 	  //putting the active nodes into the stack
   633 	  //putting the active nodes into the stack
   627 	  int lev=level[v];
   634 	  int lev=level[v];
   628 	  if ( exc > 0 && lev < n && v != t ) active[lev].push(v);
   635 	  if ( exc > 0 && lev < n && v != t ) active[lev].push(v);
   629 	}
   636 	}
   630 	break;
   637 	break;
   631       }
   638       }
   632     case GEN_FLOW:
   639     case GEN_FLOW:
   633       {
   640       {
   634 	//Counting the excess of t
   641 	//computing the excess of t
   635 	Num exc=0;
   642 	Num exc=0;
   636 	  
   643 
   637 	InEdgeIt e;
   644 	InEdgeIt e;
   638 	for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
   645 	for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
   639 	OutEdgeIt f;
   646 	OutEdgeIt f;
   640 	for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
   647 	for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
   641 	  
   648 
   642 	excess.set(t,exc);	
   649 	excess.set(t,exc);
   643 	  
   650 
   644 	break;
   651 	break;
   645       }
   652       }
   646     default:
   653     default:;
   647       break;
   654     }
   648     }
   655 
   649       
   656     preflowPreproc(fe, active, level_list, left, right);
   650     preflowPreproc( fe, active, level_list, left, right );
   657     //End of preprocessing
   651     //End of preprocessing 
   658 
   652       
   659 
   653       
       
   654     //Push/relabel on the highest level active nodes.
   660     //Push/relabel on the highest level active nodes.
   655     while ( true ) {
   661     while ( true ) {
   656       if ( b == 0 ) {
   662       if ( b == 0 ) {
   657 	if ( !what_heur && !end && k > 0 ) {
   663 	if ( !what_heur && !end && k > 0 ) {
   658 	  b=k;
   664 	  b=k;
   659 	  end=true;
   665 	  end=true;
   660 	} else break;
   666 	} else break;
   661       }
   667       }
   662 	
   668 
   663       if ( active[b].empty() ) --b; 
   669       if ( active[b].empty() ) --b;
   664       else {
   670       else {
   665 	end=false;  
   671 	end=false;
   666 	Node w=active[b].top();
   672 	Node w=active[b].top();
   667 	active[b].pop();
   673 	active[b].pop();
   668 	int newlevel=push(w,active);
   674 	int newlevel=push(w,active);
   669 	if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list, 
   675 	if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list,
   670 				     left, right, b, k, what_heur);
   676 				     left, right, b, k, what_heur);
   671 	  
   677 
   672 	++numrelabel; 
   678 	++numrelabel;
   673 	if ( numrelabel >= heur ) {
   679 	if ( numrelabel >= heur ) {
   674 	  numrelabel=0;
   680 	  numrelabel=0;
   675 	  if ( what_heur ) {
   681 	  if ( what_heur ) {
   676 	    what_heur=0;
   682 	    what_heur=0;
   677 	    heur=heur0;
   683 	    heur=heur0;
   678 	    end=false;
   684 	    end=false;
   679 	  } else {
   685 	  } else {
   680 	    what_heur=1;
   686 	    what_heur=1;
   681 	    heur=heur1;
   687 	    heur=heur1;
   682 	    b=k; 
   688 	    b=k;
   683 	  }
   689 	  }
   684 	}
   690 	}
   685       } 
   691       }
   686     } 
   692     }
   687   }
   693   }
   688 
   694 
   689 
   695 
   690 
   696 
   691   template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   697   template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   692   void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase1() 
   698   void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase1()
   693   {
   699   {
   694       
   700 
   695     int k=n-2;  //bound on the highest level under n containing a node
   701     int k=n-2;  //bound on the highest level under n containing a node
   696     int b=k;    //bound on the highest level under n of an active node
   702     int b=k;    //bound on the highest level under n of an active node
   697       
   703 
   698     VecStack active(n);
   704     VecStack active(n);
   699     level.set(s,0);
   705     level.set(s,0);
   700     std::queue<Node> bfs_queue;
   706     std::queue<Node> bfs_queue;
   701     bfs_queue.push(s);
   707     bfs_queue.push(s);
   702 	    
   708 
   703     while (!bfs_queue.empty()) {
   709     while (!bfs_queue.empty()) {
   704 	
   710 
   705       Node v=bfs_queue.front();	
   711       Node v=bfs_queue.front();
   706       bfs_queue.pop();
   712       bfs_queue.pop();
   707       int l=level[v]+1;
   713       int l=level[v]+1;
   708 	      
   714 
   709       InEdgeIt e;
   715       InEdgeIt e;
   710       for(g->first(e,v); g->valid(e); g->next(e)) {
   716       for(g->first(e,v); g->valid(e); g->next(e)) {
   711 	if ( (*capacity)[e] <= (*flow)[e] ) continue;
   717 	if ( (*capacity)[e] <= (*flow)[e] ) continue;
   712 	Node u=g->tail(e);
   718 	Node u=g->tail(e);
   713 	if ( level[u] >= n ) { 
   719 	if ( level[u] >= n ) {
   714 	  bfs_queue.push(u);
   720 	  bfs_queue.push(u);
   715 	  level.set(u, l);
   721 	  level.set(u, l);
   716 	  if ( excess[u] > 0 ) active[l].push(u);
   722 	  if ( excess[u] > 0 ) active[l].push(u);
   717 	}
   723 	}
   718       }
   724       }
   719 	
   725 
   720       OutEdgeIt f;
   726       OutEdgeIt f;
   721       for(g->first(f,v); g->valid(f); g->next(f)) {
   727       for(g->first(f,v); g->valid(f); g->next(f)) {
   722 	if ( 0 >= (*flow)[f] ) continue;
   728 	if ( 0 >= (*flow)[f] ) continue;
   723 	Node u=g->head(f);
   729 	Node u=g->head(f);
   724 	if ( level[u] >= n ) { 
   730 	if ( level[u] >= n ) {
   725 	  bfs_queue.push(u);
   731 	  bfs_queue.push(u);
   726 	  level.set(u, l);
   732 	  level.set(u, l);
   727 	  if ( excess[u] > 0 ) active[l].push(u);
   733 	  if ( excess[u] > 0 ) active[l].push(u);
   728 	}
   734 	}
   729       }
   735       }
   730     }
   736     }
   731     b=n-2;
   737     b=n-2;
   732 
   738 
   733     while ( true ) {
   739     while ( true ) {
   734 	
   740 
   735       if ( b == 0 ) break;
   741       if ( b == 0 ) break;
   736 
   742 
   737       if ( active[b].empty() ) --b; 
   743       if ( active[b].empty() ) --b;
   738       else {
   744       else {
   739 	Node w=active[b].top();
   745 	Node w=active[b].top();
   740 	active[b].pop();
   746 	active[b].pop();
   741 	int newlevel=push(w,active);	  
   747 	int newlevel=push(w,active);
   742 
   748 
   743 	//relabel
   749 	//relabel
   744 	if ( excess[w] > 0 ) {
   750 	if ( excess[w] > 0 ) {
   745 	  level.set(w,++newlevel);
   751 	  level.set(w,++newlevel);
   746 	  active[newlevel].push(w);
   752 	  active[newlevel].push(w);
   751   }
   757   }
   752 
   758 
   753 
   759 
   754 
   760 
   755   template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   761   template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   756   bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath() 
   762   bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath()
   757   {
   763   {
   758     ResGW res_graph(*g, *capacity, *flow);
   764     ResGW res_graph(*g, *capacity, *flow);
   759     bool _augment=false;
   765     bool _augment=false;
   760       
   766 
   761     //ReachedMap level(res_graph);
   767     //ReachedMap level(res_graph);
   762     FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
   768     FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
   763     BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
   769     BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
   764     bfs.pushAndSetReached(s);
   770     bfs.pushAndSetReached(s);
   765 	
   771 
   766     typename ResGW::template NodeMap<ResGWEdge> pred(res_graph); 
   772     typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
   767     pred.set(s, INVALID);
   773     pred.set(s, INVALID);
   768       
   774 
   769     typename ResGW::template NodeMap<Num> free(res_graph);
   775     typename ResGW::template NodeMap<Num> free(res_graph);
   770 	
   776 
   771     //searching for augmenting path
   777     //searching for augmenting path
   772     while ( !bfs.finished() ) { 
   778     while ( !bfs.finished() ) {
   773       ResGWOutEdgeIt e=bfs;
   779       ResGWOutEdgeIt e=bfs;
   774       if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
   780       if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
   775 	Node v=res_graph.tail(e);
   781 	Node v=res_graph.tail(e);
   776 	Node w=res_graph.head(e);
   782 	Node w=res_graph.head(e);
   777 	pred.set(w, e);
   783 	pred.set(w, e);
   778 	if (res_graph.valid(pred[v])) {
   784 	if (res_graph.valid(pred[v])) {
   779 	  free.set(w, std::min(free[v], res_graph.resCap(e)));
   785 	  free.set(w, std::min(free[v], res_graph.resCap(e)));
   780 	} else {
   786 	} else {
   781 	  free.set(w, res_graph.resCap(e)); 
   787 	  free.set(w, res_graph.resCap(e));
   782 	}
   788 	}
   783 	if (res_graph.head(e)==t) { _augment=true; break; }
   789 	if (res_graph.head(e)==t) { _augment=true; break; }
   784       }
   790       }
   785 	
   791 
   786       ++bfs;
   792       ++bfs;
   787     } //end of searching augmenting path
   793     } //end of searching augmenting path
   788 
   794 
   789     if (_augment) {
   795     if (_augment) {
   790       Node n=t;
   796       Node n=t;
   791       Num augment_value=free[t];
   797       Num augment_value=free[t];
   792       while (res_graph.valid(pred[n])) { 
   798       while (res_graph.valid(pred[n])) {
   793 	ResGWEdge e=pred[n];
   799 	ResGWEdge e=pred[n];
   794 	res_graph.augment(e, augment_value); 
   800 	res_graph.augment(e, augment_value);
   795 	n=res_graph.tail(e);
   801 	n=res_graph.tail(e);
   796       }
   802       }
   797     }
   803     }
   798 
   804 
   799     return _augment;
   805     return _augment;
   803 
   809 
   804 
   810 
   805 
   811 
   806 
   812 
   807 
   813 
   808 
       
   809 
       
   810   template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   814   template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   811   template<typename MutableGraph> 
   815   template<typename MutableGraph>
   812   bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow() 
   816   bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow()
   813   {      
   817   {
   814     typedef MutableGraph MG;
   818     typedef MutableGraph MG;
   815     bool _augment=false;
   819     bool _augment=false;
   816 
   820 
   817     ResGW res_graph(*g, *capacity, *flow);
   821     ResGW res_graph(*g, *capacity, *flow);
   818 
   822 
   819     //bfs for distances on the residual graph
   823     //bfs for distances on the residual graph
   820     //ReachedMap level(res_graph);
   824     //ReachedMap level(res_graph);
   821     FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
   825     FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
   822     BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
   826     BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
   823     bfs.pushAndSetReached(s);
   827     bfs.pushAndSetReached(s);
   824     typename ResGW::template NodeMap<int> 
   828     typename ResGW::template NodeMap<int>
   825       dist(res_graph); //filled up with 0's
   829       dist(res_graph); //filled up with 0's
   826 
   830 
   827     //F will contain the physical copy of the residual graph
   831     //F will contain the physical copy of the residual graph
   828     //with the set of edges which are on shortest paths
   832     //with the set of edges which are on shortest paths
   829     MG F;
   833     MG F;
   830     typename ResGW::template NodeMap<typename MG::Node> 
   834     typename ResGW::template NodeMap<typename MG::Node>
   831       res_graph_to_F(res_graph);
   835       res_graph_to_F(res_graph);
   832     {
   836     {
   833       typename ResGW::NodeIt n;
   837       typename ResGW::NodeIt n;
   834       for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
   838       for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
   835 	res_graph_to_F.set(n, F.addNode());
   839 	res_graph_to_F.set(n, F.addNode());
   839     typename MG::Node sF=res_graph_to_F[s];
   843     typename MG::Node sF=res_graph_to_F[s];
   840     typename MG::Node tF=res_graph_to_F[t];
   844     typename MG::Node tF=res_graph_to_F[t];
   841     typename MG::template EdgeMap<ResGWEdge> original_edge(F);
   845     typename MG::template EdgeMap<ResGWEdge> original_edge(F);
   842     typename MG::template EdgeMap<Num> residual_capacity(F);
   846     typename MG::template EdgeMap<Num> residual_capacity(F);
   843 
   847 
   844     while ( !bfs.finished() ) { 
   848     while ( !bfs.finished() ) {
   845       ResGWOutEdgeIt e=bfs;
   849       ResGWOutEdgeIt e=bfs;
   846       if (res_graph.valid(e)) {
   850       if (res_graph.valid(e)) {
   847 	if (bfs.isBNodeNewlyReached()) {
   851 	if (bfs.isBNodeNewlyReached()) {
   848 	  dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
   852 	  dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
   849 	  typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
   853 	  typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)],
       
   854 					res_graph_to_F[res_graph.head(e)]);
   850 	  original_edge.update();
   855 	  original_edge.update();
   851 	  original_edge.set(f, e);
   856 	  original_edge.set(f, e);
   852 	  residual_capacity.update();
   857 	  residual_capacity.update();
   853 	  residual_capacity.set(f, res_graph.resCap(e));
   858 	  residual_capacity.set(f, res_graph.resCap(e));
   854 	} else {
   859 	} else {
   855 	  if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
   860 	  if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
   856 	    typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
   861 	    typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)],
       
   862 					  res_graph_to_F[res_graph.head(e)]);
   857 	    original_edge.update();
   863 	    original_edge.update();
   858 	    original_edge.set(f, e);
   864 	    original_edge.set(f, e);
   859 	    residual_capacity.update();
   865 	    residual_capacity.update();
   860 	    residual_capacity.set(f, res_graph.resCap(e));
   866 	    residual_capacity.set(f, res_graph.resCap(e));
   861 	  }
   867 	  }
   874       pred.set(sF, INVALID);
   880       pred.set(sF, INVALID);
   875       //invalid iterators for sources
   881       //invalid iterators for sources
   876 
   882 
   877       typename MG::template NodeMap<Num> free(F);
   883       typename MG::template NodeMap<Num> free(F);
   878 
   884 
   879       dfs.pushAndSetReached(sF);      
   885       dfs.pushAndSetReached(sF);
   880       while (!dfs.finished()) {
   886       while (!dfs.finished()) {
   881 	++dfs;
   887 	++dfs;
   882 	if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
   888 	if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
   883 	  if (dfs.isBNodeNewlyReached()) {
   889 	  if (dfs.isBNodeNewlyReached()) {
   884 	    typename MG::Node v=F.aNode(dfs);
   890 	    typename MG::Node v=F.aNode(dfs);
   885 	    typename MG::Node w=F.bNode(dfs);
   891 	    typename MG::Node w=F.bNode(dfs);
   886 	    pred.set(w, dfs);
   892 	    pred.set(w, dfs);
   887 	    if (F.valid(pred[v])) {
   893 	    if (F.valid(pred[v])) {
   888 	      free.set(w, std::min(free[v], residual_capacity[dfs]));
   894 	      free.set(w, std::min(free[v], residual_capacity[dfs]));
   889 	    } else {
   895 	    } else {
   890 	      free.set(w, residual_capacity[dfs]); 
   896 	      free.set(w, residual_capacity[dfs]);
   891 	    }
   897 	    }
   892 	    if (w==tF) { 
   898 	    if (w==tF) {
   893 	      __augment=true; 
   899 	      __augment=true;
   894 	      _augment=true;
   900 	      _augment=true;
   895 	      break; 
   901 	      break;
   896 	    }
   902 	    }
   897 	      
   903 
   898 	  } else {
   904 	  } else {
   899 	    F.erase(/*typename MG::OutEdgeIt*/(dfs));
   905 	    F.erase(/*typename MG::OutEdgeIt*/(dfs));
   900 	  }
   906 	  }
   901 	} 
   907 	}
   902       }
   908       }
   903 
   909 
   904       if (__augment) {
   910       if (__augment) {
   905 	typename MG::Node n=tF;
   911 	typename MG::Node n=tF;
   906 	Num augment_value=free[tF];
   912 	Num augment_value=free[tF];
   907 	while (F.valid(pred[n])) { 
   913 	while (F.valid(pred[n])) {
   908 	  typename MG::Edge e=pred[n];
   914 	  typename MG::Edge e=pred[n];
   909 	  res_graph.augment(original_edge[e], augment_value); 
   915 	  res_graph.augment(original_edge[e], augment_value);
   910 	  n=F.tail(e);
   916 	  n=F.tail(e);
   911 	  if (residual_capacity[e]==augment_value) 
   917 	  if (residual_capacity[e]==augment_value)
   912 	    F.erase(e); 
   918 	    F.erase(e);
   913 	  else 
   919 	  else
   914 	    residual_capacity.set(e, residual_capacity[e]-augment_value);
   920 	    residual_capacity.set(e, residual_capacity[e]-augment_value);
   915 	}
   921 	}
   916       }
   922       }
   917 	
   923 
   918     }
   924     }
   919             
   925 
   920     return _augment;
   926     return _augment;
   921   }
   927   }
   922 
   928 
   923 
   929 
   924 
   930 
   925 
   931 
   926 
       
   927 
       
   928   template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   932   template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   929   bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2() 
   933   bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2()
   930   {
   934   {
   931     bool _augment=false;
   935     bool _augment=false;
   932 
   936 
   933     ResGW res_graph(*g, *capacity, *flow);
   937     ResGW res_graph(*g, *capacity, *flow);
   934       
   938 
   935     //ReachedMap level(res_graph);
   939     //ReachedMap level(res_graph);
   936     FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
   940     FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
   937     BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
   941     BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
   938 
   942 
   939     bfs.pushAndSetReached(s);
   943     bfs.pushAndSetReached(s);
   940     DistanceMap<ResGW> dist(res_graph);
   944     DistanceMap<ResGW> dist(res_graph);
   941     while ( !bfs.finished() ) { 
   945     while ( !bfs.finished() ) {
   942       ResGWOutEdgeIt e=bfs;
   946       ResGWOutEdgeIt e=bfs;
   943       if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
   947       if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
   944 	dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
   948 	dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
   945       }
   949       }
   946       ++bfs;
   950       ++bfs;
   947     } //computing distances from s in the residual graph
   951     } //computing distances from s in the residual graph
   948 
   952 
   949       //Subgraph containing the edges on some shortest paths
   953       //Subgraph containing the edges on some shortest paths
   950     ConstMap<typename ResGW::Node, bool> true_map(true);
   954     ConstMap<typename ResGW::Node, bool> true_map(true);
   951     typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>, 
   955     typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>,
   952       DistanceMap<ResGW> > FilterResGW;
   956       DistanceMap<ResGW> > FilterResGW;
   953     FilterResGW filter_res_graph(res_graph, true_map, dist);
   957     FilterResGW filter_res_graph(res_graph, true_map, dist);
   954 
   958 
   955     //Subgraph, which is able to delete edges which are already 
   959     //Subgraph, which is able to delete edges which are already
   956     //met by the dfs
   960     //met by the dfs
   957     typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt> 
   961     typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt>
   958       first_out_edges(filter_res_graph);
   962       first_out_edges(filter_res_graph);
   959     typename FilterResGW::NodeIt v;
   963     typename FilterResGW::NodeIt v;
   960     for(filter_res_graph.first(v); filter_res_graph.valid(v); 
   964     for(filter_res_graph.first(v); filter_res_graph.valid(v);
   961 	filter_res_graph.next(v)) 
   965 	filter_res_graph.next(v))
   962       {
   966       {
   963  	typename FilterResGW::OutEdgeIt e;
   967  	typename FilterResGW::OutEdgeIt e;
   964  	filter_res_graph.first(e, v);
   968  	filter_res_graph.first(e, v);
   965  	first_out_edges.set(v, e);
   969  	first_out_edges.set(v, e);
   966       }
   970       }
   972 
   976 
   973     while (__augment) {
   977     while (__augment) {
   974 
   978 
   975       __augment=false;
   979       __augment=false;
   976       //computing blocking flow with dfs
   980       //computing blocking flow with dfs
   977       DfsIterator< ErasingResGW, 
   981       DfsIterator< ErasingResGW,
   978 	typename ErasingResGW::template NodeMap<bool> > 
   982 	typename ErasingResGW::template NodeMap<bool> >
   979 	dfs(erasing_res_graph);
   983 	dfs(erasing_res_graph);
   980       typename ErasingResGW::
   984       typename ErasingResGW::
   981 	template NodeMap<typename ErasingResGW::OutEdgeIt> 
   985 	template NodeMap<typename ErasingResGW::OutEdgeIt>
   982 	pred(erasing_res_graph); 
   986 	pred(erasing_res_graph);
   983       pred.set(s, INVALID);
   987       pred.set(s, INVALID);
   984       //invalid iterators for sources
   988       //invalid iterators for sources
   985 
   989 
   986       typename ErasingResGW::template NodeMap<Num> 
   990       typename ErasingResGW::template NodeMap<Num>
   987 	free1(erasing_res_graph);
   991 	free1(erasing_res_graph);
   988 
   992 
   989       dfs.pushAndSetReached(
   993       dfs.pushAndSetReached
   990 			    typename ErasingResGW::Node(
   994 	///\bug hugo 0.2
   991 							typename FilterResGW::Node(
   995 	(typename ErasingResGW::Node
   992 										   typename ResGW::Node(s)
   996 	 (typename FilterResGW::Node
   993 										   )
   997 	  (typename ResGW::Node(s)
   994 							)
   998 	   )
   995 			    );
   999 	  )
       
  1000 	 );
   996       while (!dfs.finished()) {
  1001       while (!dfs.finished()) {
   997 	++dfs;
  1002 	++dfs;
   998 	if (erasing_res_graph.valid(
  1003 	if (erasing_res_graph.valid(typename ErasingResGW::OutEdgeIt(dfs)))
   999 				    typename ErasingResGW::OutEdgeIt(dfs))) 
  1004  	  {
  1000  	  { 
       
  1001   	    if (dfs.isBNodeNewlyReached()) {
  1005   	    if (dfs.isBNodeNewlyReached()) {
  1002 	  
  1006 
  1003  	      typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs);
  1007  	      typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs);
  1004  	      typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs);
  1008  	      typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs);
  1005 
  1009 
  1006  	      pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs));
  1010  	      pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs));
  1007  	      if (erasing_res_graph.valid(pred[v])) {
  1011  	      if (erasing_res_graph.valid(pred[v])) {
  1008  		free1.set(w, std::min(free1[v], res_graph.resCap(
  1012  		free1.set
  1009 								 typename ErasingResGW::OutEdgeIt(dfs))));
  1013 		  (w, std::min(free1[v], res_graph.resCap
       
  1014 			       (typename ErasingResGW::OutEdgeIt(dfs))));
  1010  	      } else {
  1015  	      } else {
  1011  		free1.set(w, res_graph.resCap(
  1016  		free1.set
  1012 					      typename ErasingResGW::OutEdgeIt(dfs))); 
  1017 		  (w, res_graph.resCap
       
  1018 		   (typename ErasingResGW::OutEdgeIt(dfs)));
  1013  	      }
  1019  	      }
  1014 	      
  1020 
  1015  	      if (w==t) { 
  1021  	      if (w==t) {
  1016  		__augment=true; 
  1022  		__augment=true;
  1017  		_augment=true;
  1023  		_augment=true;
  1018  		break; 
  1024  		break;
  1019  	      }
  1025  	      }
  1020  	    } else {
  1026  	    } else {
  1021  	      erasing_res_graph.erase(dfs);
  1027  	      erasing_res_graph.erase(dfs);
  1022 	    }
  1028 	    }
  1023 	  }
  1029 	  }
  1024       }	
  1030       }
  1025 
  1031 
  1026       if (__augment) {
  1032       if (__augment) {
  1027 	typename ErasingResGW::Node n=typename FilterResGW::Node(typename ResGW::Node(t));
  1033 	typename ErasingResGW::Node
       
  1034 	  n=typename FilterResGW::Node(typename ResGW::Node(t));
  1028 	// 	  typename ResGW::NodeMap<Num> a(res_graph);
  1035 	// 	  typename ResGW::NodeMap<Num> a(res_graph);
  1029 	// 	  typename ResGW::Node b;
  1036 	// 	  typename ResGW::Node b;
  1030 	// 	  Num j=a[b];
  1037 	// 	  Num j=a[b];
  1031 	// 	  typename FilterResGW::NodeMap<Num> a1(filter_res_graph);
  1038 	// 	  typename FilterResGW::NodeMap<Num> a1(filter_res_graph);
  1032 	// 	  typename FilterResGW::Node b1;
  1039 	// 	  typename FilterResGW::Node b1;
  1033 	// 	  Num j1=a1[b1];
  1040 	// 	  Num j1=a1[b1];
  1034 	// 	  typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);
  1041 	// 	  typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);
  1035 	// 	  typename ErasingResGW::Node b2;
  1042 	// 	  typename ErasingResGW::Node b2;
  1036 	// 	  Num j2=a2[b2];
  1043 	// 	  Num j2=a2[b2];
  1037 	Num augment_value=free1[n];
  1044 	Num augment_value=free1[n];
  1038 	while (erasing_res_graph.valid(pred[n])) { 
  1045 	while (erasing_res_graph.valid(pred[n])) {
  1039 	  typename ErasingResGW::OutEdgeIt e=pred[n];
  1046 	  typename ErasingResGW::OutEdgeIt e=pred[n];
  1040 	  res_graph.augment(e, augment_value);
  1047 	  res_graph.augment(e, augment_value);
  1041 	  n=erasing_res_graph.tail(e);
  1048 	  n=erasing_res_graph.tail(e);
  1042 	  if (res_graph.resCap(e)==0)
  1049 	  if (res_graph.resCap(e)==0)
  1043 	    erasing_res_graph.erase(e);
  1050 	    erasing_res_graph.erase(e);
  1044 	}
  1051 	}
  1045       }
  1052       }
  1046       
  1053 
  1047     } //while (__augment) 
  1054     } //while (__augment)
  1048             
  1055 
  1049     return _augment;
  1056     return _augment;
  1050   }
  1057   }
  1051 
  1058 
  1052 
  1059 
  1053 
       
  1054 
       
  1055 } //namespace hugo
  1060 } //namespace hugo
  1056 
  1061 
  1057 #endif //HUGO_MAX_FLOW_H
  1062 #endif //HUGO_MAX_FLOW_H
  1058 
  1063 
  1059 
  1064