src/lemon/preflow.h
changeset 1222 a3fb216a267d
parent 1164 80bb73097736
child 1227 01f668e3e168
equal deleted inserted replaced
4:7e427a6cdf9f 5:1dccbc6f5745
    35 
    35 
    36   ///%Preflow algorithms class.
    36   ///%Preflow algorithms class.
    37 
    37 
    38   ///This class provides an implementation of the \e preflow \e
    38   ///This class provides an implementation of the \e preflow \e
    39   ///algorithm producing a flow of maximum value in a directed
    39   ///algorithm producing a flow of maximum value in a directed
    40   ///graph. The preflow algorithms are the fastest max flow algorithms
    40   ///graph. The preflow algorithms are the fastest known max flow algorithms
    41   ///up to now. The \e source node, the \e target node, the \e
    41   ///up to now. The \e source node, the \e target node, the \e
    42   ///capacity of the edges and the \e starting \e flow value of the
    42   ///capacity of the edges and the \e starting \e flow value of the
    43   ///edges should be passed to the algorithm through the
    43   ///edges should be passed to the algorithm through the
    44   ///constructor. It is possible to change these quantities using the
    44   ///constructor. It is possible to change these quantities using the
    45   ///functions \ref setSource, \ref setTarget, \ref setCap and \ref
    45   ///functions \ref source, \ref target, \ref setCap and \ref
    46   ///setFlow.
    46   ///setFlow.
    47   ///
    47   ///
    48   ///After running \ref lemon::Preflow::phase1() "phase1()"
    48   ///After running \ref lemon::Preflow::phase1() "phase1()"
    49   ///or \ref lemon::Preflow::run() "run()", the maximal flow
    49   ///or \ref lemon::Preflow::run() "run()", the maximal flow
    50   ///value can be obtained by calling \ref flowValue(). The minimum
    50   ///value can be obtained by calling \ref flowValue(). The minimum
    53   ///the inclusionwise minimum and maximum of the minimum value cuts,
    53   ///the inclusionwise minimum and maximum of the minimum value cuts,
    54   ///resp.)
    54   ///resp.)
    55   ///
    55   ///
    56   ///\param Graph The directed graph type the algorithm runs on.
    56   ///\param Graph The directed graph type the algorithm runs on.
    57   ///\param Num The number type of the capacities and the flow values.
    57   ///\param Num The number type of the capacities and the flow values.
    58   ///\param CapMap The capacity map type.
    58   ///\param CapacityMap The capacity map type.
    59   ///\param FlowMap The flow map type.
    59   ///\param FlowMap The flow map type.
    60   ///
    60   ///
    61   ///\author Jacint Szabo 
    61   ///\author Jacint Szabo 
    62   template <typename Graph, typename Num,
    62   template <typename Graph, typename Num,
    63 	    typename CapMap=typename Graph::template EdgeMap<Num>,
    63 	    typename CapacityMap=typename Graph::template EdgeMap<Num>,
    64             typename FlowMap=typename Graph::template EdgeMap<Num> >
    64             typename FlowMap=typename Graph::template EdgeMap<Num> >
    65   class Preflow {
    65   class Preflow {
    66   protected:
    66   protected:
    67     typedef typename Graph::Node Node;
    67     typedef typename Graph::Node Node;
    68     typedef typename Graph::NodeIt NodeIt;
    68     typedef typename Graph::NodeIt NodeIt;
    71     typedef typename Graph::InEdgeIt InEdgeIt;
    71     typedef typename Graph::InEdgeIt InEdgeIt;
    72 
    72 
    73     typedef typename Graph::template NodeMap<Node> NNMap;
    73     typedef typename Graph::template NodeMap<Node> NNMap;
    74     typedef typename std::vector<Node> VecNode;
    74     typedef typename std::vector<Node> VecNode;
    75 
    75 
    76     const Graph* g;
    76     const Graph* _g;
    77     Node s;
    77     Node _source;
    78     Node t;
    78     Node _target;
    79     const CapMap* capacity;
    79     const CapacityMap* _capacity;
    80     FlowMap* flow;
    80     FlowMap* _flow;
    81     int n;      //the number of nodes of G
    81     int _node_num;      //the number of nodes of G
    82     
    82     
    83     typename Graph::template NodeMap<int> level;  
    83     typename Graph::template NodeMap<int> level;  
    84     typename Graph::template NodeMap<Num> excess;
    84     typename Graph::template NodeMap<Num> excess;
    85 
    85 
    86     // constants used for heuristics
    86     // constants used for heuristics
    89 
    89 
    90     public:
    90     public:
    91 
    91 
    92     ///Indicates the property of the starting flow map.
    92     ///Indicates the property of the starting flow map.
    93 
    93 
    94     ///Indicates the property of the starting flow map. The meanings are as follows:
    94     ///Indicates the property of the starting flow map.
       
    95     ///The meanings are as follows:
    95     ///- \c ZERO_FLOW: constant zero flow
    96     ///- \c ZERO_FLOW: constant zero flow
    96     ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to
    97     ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to
    97     ///the sum of the out-flows in every node except the \e source and
    98     ///the sum of the out-flows in every node except the \e source and
    98     ///the \e target.
    99     ///the \e target.
    99     ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at 
   100     ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at 
   109       PRE_FLOW
   110       PRE_FLOW
   110     };
   111     };
   111 
   112 
   112     ///Indicates the state of the preflow algorithm.
   113     ///Indicates the state of the preflow algorithm.
   113 
   114 
   114     ///Indicates the state of the preflow algorithm. The meanings are as follows:
   115     ///Indicates the state of the preflow algorithm.
   115     ///- \c AFTER_NOTHING: before running the algorithm or at an unspecified state.
   116     ///The meanings are as follows:
       
   117     ///- \c AFTER_NOTHING: before running the algorithm or
       
   118     ///  at an unspecified state.
   116     ///- \c AFTER_PREFLOW_PHASE_1: right after running \c phase1
   119     ///- \c AFTER_PREFLOW_PHASE_1: right after running \c phase1
   117     ///- \c AFTER_PREFLOW_PHASE_2: after running \ref phase2()
   120     ///- \c AFTER_PREFLOW_PHASE_2: after running \ref phase2()
   118     ///
   121     ///
   119     enum StatusEnum {
   122     enum StatusEnum {
   120       AFTER_NOTHING,
   123       AFTER_NOTHING,
   131 
   134 
   132     ///The constructor of the class. 
   135     ///The constructor of the class. 
   133     ///\param _G The directed graph the algorithm runs on. 
   136     ///\param _G The directed graph the algorithm runs on. 
   134     ///\param _s The source node.
   137     ///\param _s The source node.
   135     ///\param _t The target node.
   138     ///\param _t The target node.
   136     ///\param _capacity The capacity of the edges. 
   139     ///\param _cap The capacity of the edges. 
   137     ///\param _flow The flow of the edges. 
   140     ///\param _f The flow of the edges. 
   138     ///Except the graph, all of these parameters can be reset by
   141     ///Except the graph, all of these parameters can be reset by
   139     ///calling \ref setSource, \ref setTarget, \ref setCap and \ref
   142     ///calling \ref source, \ref target, \ref setCap and \ref
   140     ///setFlow, resp.
   143     ///setFlow, resp.
   141       Preflow(const Graph& _G, Node _s, Node _t, 
   144       Preflow(const Graph& _gr, Node _s, Node _t, 
   142 	      const CapMap& _capacity, FlowMap& _flow) :
   145 	      const CapacityMap& _cap, FlowMap& _f) :
   143 	g(&_G), s(_s), t(_t), capacity(&_capacity),
   146 	_g(&_gr), _source(_s), _target(_t), _capacity(&_cap),
   144 	flow(&_flow), n(countNodes(_G)), level(_G), excess(_G,0), 
   147 	_flow(&_f), _node_num(countNodes(_gr)), level(_gr), excess(_gr,0), 
   145 	flow_prop(NO_FLOW), status(AFTER_NOTHING) { }
   148 	flow_prop(NO_FLOW), status(AFTER_NOTHING) { }
   146 
   149 
   147 
   150 
   148                                                                               
   151                                                                               
   149     ///Runs the preflow algorithm.  
   152     ///Runs the preflow algorithm.  
   202     ///minCut returns a minimum cut.
   205     ///minCut returns a minimum cut.
   203     ///\warning \ref minCut(), \ref minMinCut() and \ref maxMinCut() do not
   206     ///\warning \ref minCut(), \ref minMinCut() and \ref maxMinCut() do not
   204     ///give minimum value cuts unless calling \ref phase2().
   207     ///give minimum value cuts unless calling \ref phase2().
   205     void phase1()
   208     void phase1()
   206     {
   209     {
   207       int heur0=(int)(H0*n);  //time while running 'bound decrease'
   210       int heur0=(int)(H0*_node_num);  //time while running 'bound decrease'
   208       int heur1=(int)(H1*n);  //time while running 'highest label'
   211       int heur1=(int)(H1*_node_num);  //time while running 'highest label'
   209       int heur=heur1;         //starting time interval (#of relabels)
   212       int heur=heur1;         //starting time interval (#of relabels)
   210       int numrelabel=0;
   213       int numrelabel=0;
   211 
   214 
   212       bool what_heur=1;
   215       bool what_heur=1;
   213       //It is 0 in case 'bound decrease' and 1 in case 'highest label'
   216       //It is 0 in case 'bound decrease' and 1 in case 'highest label'
   214 
   217 
   215       bool end=false;
   218       bool end=false;
   216       //Needed for 'bound decrease', true means no active 
   219       //Needed for 'bound decrease', true means no active 
   217       //nodes are above bound b.
   220       //nodes are above bound b.
   218 
   221 
   219       int k=n-2;  //bound on the highest level under n containing a node
   222       int k=_node_num-2;  //bound on the highest level under n containing a node
   220       int b=k;    //bound on the highest level under n of an active node
   223       int b=k;    //bound on the highest level under n of an active node
   221 
   224 
   222       VecNode first(n, INVALID);
   225       VecNode first(_node_num, INVALID);
   223       NNMap next(*g, INVALID);
   226       NNMap next(*_g, INVALID);
   224 
   227 
   225       NNMap left(*g, INVALID);
   228       NNMap left(*_g, INVALID);
   226       NNMap right(*g, INVALID);
   229       NNMap right(*_g, INVALID);
   227       VecNode level_list(n,INVALID);
   230       VecNode level_list(_node_num,INVALID);
   228       //List of the nodes in level i<n, set to n.
   231       //List of the nodes in level i<n, set to n.
   229 
   232 
   230       preflowPreproc(first, next, level_list, left, right);
   233       preflowPreproc(first, next, level_list, left, right);
   231 
   234 
   232       //Push/relabel on the highest level active nodes.
   235       //Push/relabel on the highest level active nodes.
   269     //   2 phase
   272     //   2 phase
   270     //   gap
   273     //   gap
   271     //   list 'level_list' on the nodes on level i implemented by hand
   274     //   list 'level_list' on the nodes on level i implemented by hand
   272     //   stack 'active' on the active nodes on level i      
   275     //   stack 'active' on the active nodes on level i      
   273     //   runs heuristic 'highest label' for H1*n relabels
   276     //   runs heuristic 'highest label' for H1*n relabels
   274     //   runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
   277     //   runs heuristic 'bound decrease' for H0*n relabels,
       
   278     //        starts with 'highest label'
   275     //   Parameters H0 and H1 are initialized to 20 and 1.
   279     //   Parameters H0 and H1 are initialized to 20 and 1.
   276 
   280 
   277 
   281 
   278     ///Runs the second phase of the preflow algorithm.
   282     ///Runs the second phase of the preflow algorithm.
   279 
   283 
   285     ///maxMinCut return the inclusionwise minimum and maximum cuts of
   289     ///maxMinCut return the inclusionwise minimum and maximum cuts of
   286     ///minimum value, resp.  \pre \ref phase1 must be called before.
   290     ///minimum value, resp.  \pre \ref phase1 must be called before.
   287     void phase2()
   291     void phase2()
   288     {
   292     {
   289 
   293 
   290       int k=n-2;  //bound on the highest level under n containing a node
   294       int k=_node_num-2;  //bound on the highest level under n containing a node
   291       int b=k;    //bound on the highest level under n of an active node
   295       int b=k;    //bound on the highest level under n of an active node
   292 
   296 
   293     
   297     
   294       VecNode first(n, INVALID);
   298       VecNode first(_node_num, INVALID);
   295       NNMap next(*g, INVALID); 
   299       NNMap next(*_g, INVALID); 
   296       level.set(s,0);
   300       level.set(_source,0);
   297       std::queue<Node> bfs_queue;
   301       std::queue<Node> bfs_queue;
   298       bfs_queue.push(s);
   302       bfs_queue.push(_source);
   299 
   303 
   300       while ( !bfs_queue.empty() ) {
   304       while ( !bfs_queue.empty() ) {
   301 
   305 
   302 	Node v=bfs_queue.front();
   306 	Node v=bfs_queue.front();
   303 	bfs_queue.pop();
   307 	bfs_queue.pop();
   304 	int l=level[v]+1;
   308 	int l=level[v]+1;
   305 
   309 
   306 	for(InEdgeIt e(*g,v); e!=INVALID; ++e) {
   310 	for(InEdgeIt e(*_g,v); e!=INVALID; ++e) {
   307 	  if ( (*capacity)[e] <= (*flow)[e] ) continue;
   311 	  if ( (*_capacity)[e] <= (*_flow)[e] ) continue;
   308 	  Node u=g->source(e);
   312 	  Node u=_g->source(e);
   309 	  if ( level[u] >= n ) {
   313 	  if ( level[u] >= _node_num ) {
   310 	    bfs_queue.push(u);
   314 	    bfs_queue.push(u);
   311 	    level.set(u, l);
   315 	    level.set(u, l);
   312 	    if ( excess[u] > 0 ) {
   316 	    if ( excess[u] > 0 ) {
   313 	      next.set(u,first[l]);
   317 	      next.set(u,first[l]);
   314 	      first[l]=u;
   318 	      first[l]=u;
   315 	    }
   319 	    }
   316 	  }
   320 	  }
   317 	}
   321 	}
   318 
   322 
   319 	for(OutEdgeIt e(*g,v); e!=INVALID; ++e) {
   323 	for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) {
   320 	  if ( 0 >= (*flow)[e] ) continue;
   324 	  if ( 0 >= (*_flow)[e] ) continue;
   321 	  Node u=g->target(e);
   325 	  Node u=_g->target(e);
   322 	  if ( level[u] >= n ) {
   326 	  if ( level[u] >= _node_num ) {
   323 	    bfs_queue.push(u);
   327 	    bfs_queue.push(u);
   324 	    level.set(u, l);
   328 	    level.set(u, l);
   325 	    if ( excess[u] > 0 ) {
   329 	    if ( excess[u] > 0 ) {
   326 	      next.set(u,first[l]);
   330 	      next.set(u,first[l]);
   327 	      first[l]=u;
   331 	      first[l]=u;
   328 	    }
   332 	    }
   329 	  }
   333 	  }
   330 	}
   334 	}
   331       }
   335       }
   332       b=n-2;
   336       b=_node_num-2;
   333 
   337 
   334       while ( true ) {
   338       while ( true ) {
   335 
   339 
   336 	if ( b == 0 ) break;
   340 	if ( b == 0 ) break;
   337 	if ( first[b]==INVALID ) --b;
   341 	if ( first[b]==INVALID ) --b;
   357 
   361 
   358     /// Returns the value of the maximum flow by returning the excess
   362     /// Returns the value of the maximum flow by returning the excess
   359     /// of the target node \c t. This value equals to the value of
   363     /// of the target node \c t. This value equals to the value of
   360     /// the maximum flow already after running \ref phase1.
   364     /// the maximum flow already after running \ref phase1.
   361     Num flowValue() const {
   365     Num flowValue() const {
   362       return excess[t];
   366       return excess[_target];
   363     }
   367     }
   364 
   368 
   365 
   369 
   366     ///Returns a minimum value cut.
   370     ///Returns a minimum value cut.
   367 
   371 
   373     ///be initialized to false.
   377     ///be initialized to false.
   374     template<typename _CutMap>
   378     template<typename _CutMap>
   375     void minCut(_CutMap& M) const {
   379     void minCut(_CutMap& M) const {
   376       switch ( status ) {
   380       switch ( status ) {
   377 	case AFTER_PREFLOW_PHASE_1:
   381 	case AFTER_PREFLOW_PHASE_1:
   378 	for(NodeIt v(*g); v!=INVALID; ++v) {
   382 	for(NodeIt v(*_g); v!=INVALID; ++v) {
   379 	  if (level[v] < n) {
   383 	  if (level[v] < _node_num) {
   380 	    M.set(v, false);
   384 	    M.set(v, false);
   381 	  } else {
   385 	  } else {
   382 	    M.set(v, true);
   386 	    M.set(v, true);
   383 	  }
   387 	  }
   384 	}
   388 	}
   400     ///phase2 should already be run.
   404     ///phase2 should already be run.
   401     template<typename _CutMap>
   405     template<typename _CutMap>
   402     void minMinCut(_CutMap& M) const {
   406     void minMinCut(_CutMap& M) const {
   403 
   407 
   404       std::queue<Node> queue;
   408       std::queue<Node> queue;
   405       M.set(s,true);
   409       M.set(_source,true);
   406       queue.push(s);
   410       queue.push(s);
   407       
   411       
   408       while (!queue.empty()) {
   412       while (!queue.empty()) {
   409 	Node w=queue.front();
   413 	Node w=queue.front();
   410 	queue.pop();
   414 	queue.pop();
   411 	
   415 	
   412 	for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
   416 	for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
   413 	  Node v=g->target(e);
   417 	  Node v=_g->target(e);
   414 	  if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
   418 	  if (!M[v] && (*_flow)[e] < (*_capacity)[e] ) {
   415 	    queue.push(v);
   419 	    queue.push(v);
   416 	    M.set(v, true);
   420 	    M.set(v, true);
   417 	  }
   421 	  }
   418 	}
   422 	}
   419 	
   423 	
   420 	for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
   424 	for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
   421 	  Node v=g->source(e);
   425 	  Node v=_g->source(e);
   422 	  if (!M[v] && (*flow)[e] > 0 ) {
   426 	  if (!M[v] && (*_flow)[e] > 0 ) {
   423 	    queue.push(v);
   427 	    queue.push(v);
   424 	    M.set(v, true);
   428 	    M.set(v, true);
   425 	  }
   429 	  }
   426 	}
   430 	}
   427       }
   431       }
   434     ///backward bfs from the target node \c t in the residual graph.
   438     ///backward bfs from the target node \c t in the residual graph.
   435     ///\pre \ref phase2() or run() should already be run.
   439     ///\pre \ref phase2() or run() should already be run.
   436     template<typename _CutMap>
   440     template<typename _CutMap>
   437     void maxMinCut(_CutMap& M) const {
   441     void maxMinCut(_CutMap& M) const {
   438 
   442 
   439       for(NodeIt v(*g) ; v!=INVALID; ++v) M.set(v, true);
   443       for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true);
   440 
   444 
   441       std::queue<Node> queue;
   445       std::queue<Node> queue;
   442 
   446 
   443       M.set(t,false);
   447       M.set(_target,false);
   444       queue.push(t);
   448       queue.push(_target);
   445 
   449 
   446       while (!queue.empty()) {
   450       while (!queue.empty()) {
   447         Node w=queue.front();
   451         Node w=queue.front();
   448 	queue.pop();
   452 	queue.pop();
   449 
   453 
   450 	for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
   454 	for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
   451 	  Node v=g->source(e);
   455 	  Node v=_g->source(e);
   452 	  if (M[v] && (*flow)[e] < (*capacity)[e] ) {
   456 	  if (M[v] && (*_flow)[e] < (*_capacity)[e] ) {
   453 	    queue.push(v);
   457 	    queue.push(v);
   454 	    M.set(v, false);
   458 	    M.set(v, false);
   455 	  }
   459 	  }
   456 	}
   460 	}
   457 
   461 
   458 	for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
   462 	for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
   459 	  Node v=g->target(e);
   463 	  Node v=_g->target(e);
   460 	  if (M[v] && (*flow)[e] > 0 ) {
   464 	  if (M[v] && (*_flow)[e] > 0 ) {
   461 	    queue.push(v);
   465 	    queue.push(v);
   462 	    M.set(v, false);
   466 	    M.set(v, false);
   463 	  }
   467 	  }
   464 	}
   468 	}
   465       }
   469       }
   467 
   471 
   468     ///Sets the source node to \c _s.
   472     ///Sets the source node to \c _s.
   469 
   473 
   470     ///Sets the source node to \c _s.
   474     ///Sets the source node to \c _s.
   471     /// 
   475     /// 
   472     void setSource(Node _s) { 
   476     void source(Node _s) { 
   473       s=_s; 
   477       _source=_s; 
   474       if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
   478       if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
   475       status=AFTER_NOTHING; 
   479       status=AFTER_NOTHING; 
   476     }
   480     }
   477 
   481 
       
   482     ///Returns the source node.
       
   483 
       
   484     ///Returns the source node.
       
   485     /// 
       
   486     Node source() const { 
       
   487       return _source;
       
   488     }
       
   489 
   478     ///Sets the target node to \c _t.
   490     ///Sets the target node to \c _t.
   479 
   491 
   480     ///Sets the target node to \c _t.
   492     ///Sets the target node to \c _t.
   481     ///
   493     ///
   482     void setTarget(Node _t) { 
   494     void target(Node _t) { 
   483       t=_t; 
   495       _target=_t; 
   484       if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
   496       if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
   485       status=AFTER_NOTHING; 
   497       status=AFTER_NOTHING; 
   486     }
   498     }
   487 
   499 
       
   500     ///Returns the target node.
       
   501 
       
   502     ///Returns the target node.
       
   503     /// 
       
   504     Node target() const { 
       
   505       return _target;
       
   506     }
       
   507 
   488     /// Sets the edge map of the capacities to _cap.
   508     /// Sets the edge map of the capacities to _cap.
   489 
   509 
   490     /// Sets the edge map of the capacities to _cap.
   510     /// Sets the edge map of the capacities to _cap.
   491     /// 
   511     /// 
   492     void setCap(const CapMap& _cap) { 
   512     void capacityMap(const CapacityMap& _cap) { 
   493       capacity=&_cap; 
   513       _capacity=&_cap; 
   494       status=AFTER_NOTHING; 
   514       status=AFTER_NOTHING; 
       
   515     }
       
   516     /// Returns a reference to to capacity map.
       
   517 
       
   518     /// Returns a reference to to capacity map.
       
   519     /// 
       
   520     const CapacityMap &capacityMap() const { 
       
   521       return *_capacity;
   495     }
   522     }
   496 
   523 
   497     /// Sets the edge map of the flows to _flow.
   524     /// Sets the edge map of the flows to _flow.
   498 
   525 
   499     /// Sets the edge map of the flows to _flow.
   526     /// Sets the edge map of the flows to _flow.
   500     /// 
   527     /// 
   501     void setFlow(FlowMap& _flow) { 
   528     void flowMap(FlowMap& _f) { 
   502       flow=&_flow; 
   529       _flow=&_f; 
   503       flow_prop=NO_FLOW;
   530       flow_prop=NO_FLOW;
   504       status=AFTER_NOTHING; 
   531       status=AFTER_NOTHING; 
   505     }
   532     }
   506 
   533      
       
   534     /// Returns a reference to to flow map.
       
   535 
       
   536     /// Returns a reference to to flow map.
       
   537     /// 
       
   538     const FlowMap &flowMap() const { 
       
   539       return *_flow;
       
   540     }
   507 
   541 
   508   private:
   542   private:
   509 
   543 
   510     int push(Node w, NNMap& next, VecNode& first) {
   544     int push(Node w, NNMap& next, VecNode& first) {
   511 
   545 
   512       int lev=level[w];
   546       int lev=level[w];
   513       Num exc=excess[w];
   547       Num exc=excess[w];
   514       int newlevel=n;       //bound on the next level of w
   548       int newlevel=_node_num;       //bound on the next level of w
   515 
   549 
   516       for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
   550       for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
   517 	if ( (*flow)[e] >= (*capacity)[e] ) continue;
   551 	if ( (*_flow)[e] >= (*_capacity)[e] ) continue;
   518 	Node v=g->target(e);
   552 	Node v=_g->target(e);
   519 
   553 
   520 	if( lev > level[v] ) { //Push is allowed now
   554 	if( lev > level[v] ) { //Push is allowed now
   521 	  
   555 	  
   522 	  if ( excess[v]<=0 && v!=t && v!=s ) {
   556 	  if ( excess[v]<=0 && v!=_target && v!=_source ) {
   523 	    next.set(v,first[level[v]]);
   557 	    next.set(v,first[level[v]]);
   524 	    first[level[v]]=v;
   558 	    first[level[v]]=v;
   525 	  }
   559 	  }
   526 
   560 
   527 	  Num cap=(*capacity)[e];
   561 	  Num cap=(*_capacity)[e];
   528 	  Num flo=(*flow)[e];
   562 	  Num flo=(*_flow)[e];
   529 	  Num remcap=cap-flo;
   563 	  Num remcap=cap-flo;
   530 	  
   564 	  
   531 	  if ( remcap >= exc ) { //A nonsaturating push.
   565 	  if ( remcap >= exc ) { //A nonsaturating push.
   532 	    
   566 	    
   533 	    flow->set(e, flo+exc);
   567 	    _flow->set(e, flo+exc);
   534 	    excess.set(v, excess[v]+exc);
   568 	    excess.set(v, excess[v]+exc);
   535 	    exc=0;
   569 	    exc=0;
   536 	    break;
   570 	    break;
   537 
   571 
   538 	  } else { //A saturating push.
   572 	  } else { //A saturating push.
   539 	    flow->set(e, cap);
   573 	    _flow->set(e, cap);
   540 	    excess.set(v, excess[v]+remcap);
   574 	    excess.set(v, excess[v]+remcap);
   541 	    exc-=remcap;
   575 	    exc-=remcap;
   542 	  }
   576 	  }
   543 	} else if ( newlevel > level[v] ) newlevel = level[v];
   577 	} else if ( newlevel > level[v] ) newlevel = level[v];
   544       } //for out edges wv
   578       } //for out edges wv
   545 
   579 
   546       if ( exc > 0 ) {
   580       if ( exc > 0 ) {
   547 	for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
   581 	for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
   548 	  
   582 	  
   549 	  if( (*flow)[e] <= 0 ) continue;
   583 	  if( (*_flow)[e] <= 0 ) continue;
   550 	  Node v=g->source(e);
   584 	  Node v=_g->source(e);
   551 
   585 
   552 	  if( lev > level[v] ) { //Push is allowed now
   586 	  if( lev > level[v] ) { //Push is allowed now
   553 
   587 
   554 	    if ( excess[v]<=0 && v!=t && v!=s ) {
   588 	    if ( excess[v]<=0 && v!=_target && v!=_source ) {
   555 	      next.set(v,first[level[v]]);
   589 	      next.set(v,first[level[v]]);
   556 	      first[level[v]]=v;
   590 	      first[level[v]]=v;
   557 	    }
   591 	    }
   558 
   592 
   559 	    Num flo=(*flow)[e];
   593 	    Num flo=(*_flow)[e];
   560 
   594 
   561 	    if ( flo >= exc ) { //A nonsaturating push.
   595 	    if ( flo >= exc ) { //A nonsaturating push.
   562 
   596 
   563 	      flow->set(e, flo-exc);
   597 	      _flow->set(e, flo-exc);
   564 	      excess.set(v, excess[v]+exc);
   598 	      excess.set(v, excess[v]+exc);
   565 	      exc=0;
   599 	      exc=0;
   566 	      break;
   600 	      break;
   567 	    } else {  //A saturating push.
   601 	    } else {  //A saturating push.
   568 
   602 
   569 	      excess.set(v, excess[v]+flo);
   603 	      excess.set(v, excess[v]+flo);
   570 	      exc-=flo;
   604 	      exc-=flo;
   571 	      flow->set(e,0);
   605 	      _flow->set(e,0);
   572 	    }
   606 	    }
   573 	  } else if ( newlevel > level[v] ) newlevel = level[v];
   607 	  } else if ( newlevel > level[v] ) newlevel = level[v];
   574 	} //for in edges vw
   608 	} //for in edges vw
   575 
   609 
   576       } // if w still has excess after the out edge for cycle
   610       } // if w still has excess after the out edge for cycle
   583     
   617     
   584     
   618     
   585     void preflowPreproc(VecNode& first, NNMap& next, 
   619     void preflowPreproc(VecNode& first, NNMap& next, 
   586 			VecNode& level_list, NNMap& left, NNMap& right)
   620 			VecNode& level_list, NNMap& left, NNMap& right)
   587     {
   621     {
   588       for(NodeIt v(*g); v!=INVALID; ++v) level.set(v,n);
   622       for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num);
   589       std::queue<Node> bfs_queue;
   623       std::queue<Node> bfs_queue;
   590       
   624       
   591       if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
   625       if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
   592 	//Reverse_bfs from t in the residual graph,
   626 	//Reverse_bfs from t in the residual graph,
   593 	//to find the starting level.
   627 	//to find the starting level.
   594 	level.set(t,0);
   628 	level.set(_target,0);
   595 	bfs_queue.push(t);
   629 	bfs_queue.push(_target);
   596 	
   630 	
   597 	while ( !bfs_queue.empty() ) {
   631 	while ( !bfs_queue.empty() ) {
   598 	  
   632 	  
   599 	  Node v=bfs_queue.front();
   633 	  Node v=bfs_queue.front();
   600 	  bfs_queue.pop();
   634 	  bfs_queue.pop();
   601 	  int l=level[v]+1;
   635 	  int l=level[v]+1;
   602 	  
   636 	  
   603 	  for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) {
   637 	  for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
   604 	    if ( (*capacity)[e] <= (*flow)[e] ) continue;
   638 	    if ( (*_capacity)[e] <= (*_flow)[e] ) continue;
   605 	    Node w=g->source(e);
   639 	    Node w=_g->source(e);
   606 	    if ( level[w] == n && w != s ) {
   640 	    if ( level[w] == _node_num && w != _source ) {
   607 	      bfs_queue.push(w);
   641 	      bfs_queue.push(w);
   608 	      Node z=level_list[l];
   642 	      Node z=level_list[l];
   609 	      if ( z!=INVALID ) left.set(z,w);
   643 	      if ( z!=INVALID ) left.set(z,w);
   610 	      right.set(w,z);
   644 	      right.set(w,z);
   611 	      level_list[l]=w;
   645 	      level_list[l]=w;
   612 	      level.set(w, l);
   646 	      level.set(w, l);
   613 	    }
   647 	    }
   614 	  }
   648 	  }
   615 	  
   649 	  
   616 	  for(OutEdgeIt e(*g,v) ; e!=INVALID; ++e) {
   650 	  for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
   617 	    if ( 0 >= (*flow)[e] ) continue;
   651 	    if ( 0 >= (*_flow)[e] ) continue;
   618 	    Node w=g->target(e);
   652 	    Node w=_g->target(e);
   619 	    if ( level[w] == n && w != s ) {
   653 	    if ( level[w] == _node_num && w != _source ) {
   620 	      bfs_queue.push(w);
   654 	      bfs_queue.push(w);
   621 	      Node z=level_list[l];
   655 	      Node z=level_list[l];
   622 	      if ( z!=INVALID ) left.set(z,w);
   656 	      if ( z!=INVALID ) left.set(z,w);
   623 	      right.set(w,z);
   657 	      right.set(w,z);
   624 	      level_list[l]=w;
   658 	      level_list[l]=w;
   629       } //if
   663       } //if
   630 
   664 
   631 
   665 
   632       switch (flow_prop) {
   666       switch (flow_prop) {
   633 	case NO_FLOW:  
   667 	case NO_FLOW:  
   634 	for(EdgeIt e(*g); e!=INVALID; ++e) flow->set(e,0);
   668 	for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0);
   635 	case ZERO_FLOW:
   669 	case ZERO_FLOW:
   636 	for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0);
   670 	for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
   637 	
   671 	
   638 	//Reverse_bfs from t, to find the starting level.
   672 	//Reverse_bfs from t, to find the starting level.
   639 	level.set(t,0);
   673 	level.set(_target,0);
   640 	bfs_queue.push(t);
   674 	bfs_queue.push(_target);
   641 	
   675 	
   642 	while ( !bfs_queue.empty() ) {
   676 	while ( !bfs_queue.empty() ) {
   643 	  
   677 	  
   644 	  Node v=bfs_queue.front();
   678 	  Node v=bfs_queue.front();
   645 	  bfs_queue.pop();
   679 	  bfs_queue.pop();
   646 	  int l=level[v]+1;
   680 	  int l=level[v]+1;
   647 	  
   681 	  
   648 	  for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) {
   682 	  for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
   649 	    Node w=g->source(e);
   683 	    Node w=_g->source(e);
   650 	    if ( level[w] == n && w != s ) {
   684 	    if ( level[w] == _node_num && w != _source ) {
   651 	      bfs_queue.push(w);
   685 	      bfs_queue.push(w);
   652 	      Node z=level_list[l];
   686 	      Node z=level_list[l];
   653 	      if ( z!=INVALID ) left.set(z,w);
   687 	      if ( z!=INVALID ) left.set(z,w);
   654 	      right.set(w,z);
   688 	      right.set(w,z);
   655 	      level_list[l]=w;
   689 	      level_list[l]=w;
   657 	    }
   691 	    }
   658 	  }
   692 	  }
   659 	}
   693 	}
   660 	
   694 	
   661 	//the starting flow
   695 	//the starting flow
   662 	for(OutEdgeIt e(*g,s) ; e!=INVALID; ++e) {
   696 	for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
   663 	  Num c=(*capacity)[e];
   697 	  Num c=(*_capacity)[e];
   664 	  if ( c <= 0 ) continue;
   698 	  if ( c <= 0 ) continue;
   665 	  Node w=g->target(e);
   699 	  Node w=_g->target(e);
   666 	  if ( level[w] < n ) {
   700 	  if ( level[w] < _node_num ) {
   667 	    if ( excess[w] <= 0 && w!=t ) { //putting into the stack
   701 	    if ( excess[w] <= 0 && w!=_target ) { //putting into the stack
   668 	      next.set(w,first[level[w]]);
   702 	      next.set(w,first[level[w]]);
   669 	      first[level[w]]=w;
   703 	      first[level[w]]=w;
   670 	    }
   704 	    }
   671 	    flow->set(e, c);
   705 	    _flow->set(e, c);
   672 	    excess.set(w, excess[w]+c);
   706 	    excess.set(w, excess[w]+c);
   673 	  }
   707 	  }
   674 	}
   708 	}
   675 	break;
   709 	break;
   676 
   710 
   677 	case GEN_FLOW:
   711 	case GEN_FLOW:
   678 	for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0);
   712 	for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
   679 	{
   713 	{
   680 	  Num exc=0;
   714 	  Num exc=0;
   681 	  for(InEdgeIt e(*g,t) ; e!=INVALID; ++e) exc+=(*flow)[e];
   715 	  for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e];
   682 	  for(OutEdgeIt e(*g,t) ; e!=INVALID; ++e) exc-=(*flow)[e];
   716 	  for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e];
   683 	  excess.set(t,exc);
   717 	  excess.set(_target,exc);
   684 	}
   718 	}
   685 
   719 
   686 	//the starting flow
   720 	//the starting flow
   687 	for(OutEdgeIt e(*g,s); e!=INVALID; ++e)	{
   721 	for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e)	{
   688 	  Num rem=(*capacity)[e]-(*flow)[e];
   722 	  Num rem=(*_capacity)[e]-(*_flow)[e];
   689 	  if ( rem <= 0 ) continue;
   723 	  if ( rem <= 0 ) continue;
   690 	  Node w=g->target(e);
   724 	  Node w=_g->target(e);
   691 	  if ( level[w] < n ) {
   725 	  if ( level[w] < _node_num ) {
   692 	    if ( excess[w] <= 0 && w!=t ) { //putting into the stack
   726 	    if ( excess[w] <= 0 && w!=_target ) { //putting into the stack
   693 	      next.set(w,first[level[w]]);
   727 	      next.set(w,first[level[w]]);
   694 	      first[level[w]]=w;
   728 	      first[level[w]]=w;
   695 	    }   
   729 	    }   
   696 	    flow->set(e, (*capacity)[e]);
   730 	    _flow->set(e, (*_capacity)[e]);
   697 	    excess.set(w, excess[w]+rem);
   731 	    excess.set(w, excess[w]+rem);
   698 	  }
   732 	  }
   699 	}
   733 	}
   700 	
   734 	
   701 	for(InEdgeIt e(*g,s); e!=INVALID; ++e) {
   735 	for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) {
   702 	  if ( (*flow)[e] <= 0 ) continue;
   736 	  if ( (*_flow)[e] <= 0 ) continue;
   703 	  Node w=g->source(e);
   737 	  Node w=_g->source(e);
   704 	  if ( level[w] < n ) {
   738 	  if ( level[w] < _node_num ) {
   705 	    if ( excess[w] <= 0 && w!=t ) {
   739 	    if ( excess[w] <= 0 && w!=_target ) {
   706 	      next.set(w,first[level[w]]);
   740 	      next.set(w,first[level[w]]);
   707 	      first[level[w]]=w;
   741 	      first[level[w]]=w;
   708 	    }  
   742 	    }  
   709 	    excess.set(w, excess[w]+(*flow)[e]);
   743 	    excess.set(w, excess[w]+(*_flow)[e]);
   710 	    flow->set(e, 0);
   744 	    _flow->set(e, 0);
   711 	  }
   745 	  }
   712 	}
   746 	}
   713 	break;
   747 	break;
   714 
   748 
   715 	case PRE_FLOW:	
   749 	case PRE_FLOW:	
   716 	//the starting flow
   750 	//the starting flow
   717 	for(OutEdgeIt e(*g,s) ; e!=INVALID; ++e) {
   751 	for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
   718 	  Num rem=(*capacity)[e]-(*flow)[e];
   752 	  Num rem=(*_capacity)[e]-(*_flow)[e];
   719 	  if ( rem <= 0 ) continue;
   753 	  if ( rem <= 0 ) continue;
   720 	  Node w=g->target(e);
   754 	  Node w=_g->target(e);
   721 	  if ( level[w] < n ) flow->set(e, (*capacity)[e]);
   755 	  if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]);
   722 	}
   756 	}
   723 	
   757 	
   724 	for(InEdgeIt e(*g,s) ; e!=INVALID; ++e) {
   758 	for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
   725 	  if ( (*flow)[e] <= 0 ) continue;
   759 	  if ( (*_flow)[e] <= 0 ) continue;
   726 	  Node w=g->source(e);
   760 	  Node w=_g->source(e);
   727 	  if ( level[w] < n ) flow->set(e, 0);
   761 	  if ( level[w] < _node_num ) _flow->set(e, 0);
   728 	}
   762 	}
   729 	
   763 	
   730 	//computing the excess
   764 	//computing the excess
   731 	for(NodeIt w(*g); w!=INVALID; ++w) {
   765 	for(NodeIt w(*_g); w!=INVALID; ++w) {
   732 	  Num exc=0;
   766 	  Num exc=0;
   733 	  for(InEdgeIt e(*g,w); e!=INVALID; ++e) exc+=(*flow)[e];
   767 	  for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e];
   734 	  for(OutEdgeIt e(*g,w); e!=INVALID; ++e) exc-=(*flow)[e];
   768 	  for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e];
   735 	  excess.set(w,exc);
   769 	  excess.set(w,exc);
   736 	  
   770 	  
   737 	  //putting the active nodes into the stack
   771 	  //putting the active nodes into the stack
   738 	  int lev=level[w];
   772 	  int lev=level[w];
   739 	    if ( exc > 0 && lev < n && Node(w) != t ) {
   773 	    if ( exc > 0 && lev < _node_num && Node(w) != _target ) {
   740 	      next.set(w,first[lev]);
   774 	      next.set(w,first[lev]);
   741 	      first[lev]=w;
   775 	      first[lev]=w;
   742 	    }
   776 	    }
   743 	}
   777 	}
   744 	break;
   778 	break;
   778 
   812 
   779 	//gapping starts
   813 	//gapping starts
   780 	for (int i=lev; i!=k ; ) {
   814 	for (int i=lev; i!=k ; ) {
   781 	  Node v=level_list[++i];
   815 	  Node v=level_list[++i];
   782 	  while ( v!=INVALID ) {
   816 	  while ( v!=INVALID ) {
   783 	    level.set(v,n);
   817 	    level.set(v,_node_num);
   784 	    v=right[v];
   818 	    v=right[v];
   785 	  }
   819 	  }
   786 	  level_list[i]=INVALID;
   820 	  level_list[i]=INVALID;
   787 	  if ( !what_heur ) first[i]=INVALID;
   821 	  if ( !what_heur ) first[i]=INVALID;
   788 	}
   822 	}
   789 
   823 
   790 	level.set(w,n);
   824 	level.set(w,_node_num);
   791 	b=lev-1;
   825 	b=lev-1;
   792 	k=b;
   826 	k=b;
   793 	//gapping ends
   827 	//gapping ends
   794 
   828 
   795       } else {
   829       } else {
   796 
   830 
   797 	if ( newlevel == n ) level.set(w,n);
   831 	if ( newlevel == _node_num ) level.set(w,_node_num);
   798 	else {
   832 	else {
   799 	  level.set(w,++newlevel);
   833 	  level.set(w,++newlevel);
   800 	  next.set(w,first[newlevel]);
   834 	  next.set(w,first[newlevel]);
   801 	  first[newlevel]=w;
   835 	  first[newlevel]=w;
   802 	  if ( what_heur ) b=newlevel;
   836 	  if ( what_heur ) b=newlevel;