src/hugo/max_flow.h
changeset 753 f5382a084c07
parent 745 d976ba609099
child 757 8680351d0c28
equal deleted inserted replaced
2:b943fe5279af 3:da19b45181c9
     1 // -*- C++ -*-
     1 // -*- C++ -*-
     2 #ifndef HUGO_MAX_FLOW_NO_STACK_H
     2 #ifndef HUGO_MAX_FLOW_H
     3 #define HUGO_MAX_FLOW_NO_STACK_H
     3 #define HUGO_MAX_FLOW_H
     4 
     4 
     5 #include <vector>
     5 #include <vector>
     6 #include <queue>
     6 #include <queue>
     7 //#include <stack>
       
     8 
     7 
     9 #include <hugo/graph_wrapper.h>
     8 #include <hugo/graph_wrapper.h>
    10 #include <hugo/invalid.h>
     9 #include <hugo/invalid.h>
    11 #include <hugo/maps.h>
    10 #include <hugo/maps.h>
    12 
    11 
    13 /// \file
    12 /// \file
    14 /// \brief The same as max_flow.h, but without using stl stack for the active nodes. Only for test.
       
    15 /// \ingroup galgs
    13 /// \ingroup galgs
    16 
    14 
    17 namespace hugo {
    15 namespace hugo {
    18 
    16 
    19   /// \addtogroup galgs
    17   /// \addtogroup galgs
    49     typedef typename Graph::NodeIt NodeIt;
    47     typedef typename Graph::NodeIt NodeIt;
    50     typedef typename Graph::EdgeIt EdgeIt;
    48     typedef typename Graph::EdgeIt EdgeIt;
    51     typedef typename Graph::OutEdgeIt OutEdgeIt;
    49     typedef typename Graph::OutEdgeIt OutEdgeIt;
    52     typedef typename Graph::InEdgeIt InEdgeIt;
    50     typedef typename Graph::InEdgeIt InEdgeIt;
    53 
    51 
    54     //    typedef typename std::vector<std::stack<Node> > VecStack;
       
    55     typedef typename std::vector<Node> VecFirst;
    52     typedef typename std::vector<Node> VecFirst;
    56     typedef typename Graph::template NodeMap<Node> NNMap;
    53     typedef typename Graph::template NodeMap<Node> NNMap;
    57     typedef typename std::vector<Node> VecNode;
    54     typedef typename std::vector<Node> VecNode;
    58 
    55 
    59     const Graph* g;
    56     const Graph* g;
    64     int n;      //the number of nodes of G
    61     int n;      //the number of nodes of G
    65     typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;   
    62     typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;   
    66     //typedef ExpResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
    63     //typedef ExpResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
    67     typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
    64     typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
    68     typedef typename ResGW::Edge ResGWEdge;
    65     typedef typename ResGW::Edge ResGWEdge;
    69     //typedef typename ResGW::template NodeMap<bool> ReachedMap;
       
    70     typedef typename Graph::template NodeMap<int> ReachedMap;
    66     typedef typename Graph::template NodeMap<int> ReachedMap;
    71 
    67 
    72 
    68 
    73     //level works as a bool map in augmenting path algorithms and is
    69     //level works as a bool map in augmenting path algorithms and is
    74     //used by bfs for storing reached information.  In preflow, it
    70     //used by bfs for storing reached information.  In preflow, it
   108       AFTER_FAST_AUGMENTING, 
   104       AFTER_FAST_AUGMENTING, 
   109       AFTER_PRE_FLOW_PHASE_1,      
   105       AFTER_PRE_FLOW_PHASE_1,      
   110       AFTER_PRE_FLOW_PHASE_2
   106       AFTER_PRE_FLOW_PHASE_2
   111     };
   107     };
   112 
   108 
   113     /// Don not needle this flag only if necessary.
   109     /// Do not needle this flag only if necessary.
   114     StatusEnum status;
   110     StatusEnum status;
   115 
   111 
   116 //     int number_of_augmentations;
   112 //     int number_of_augmentations;
   117 
   113 
   118 
   114 
   186     ///Runs the first phase of the preflow algorithm.
   182     ///Runs the first phase of the preflow algorithm.
   187 
   183 
   188     ///The preflow algorithm consists of two phases, this method runs the
   184     ///The preflow algorithm consists of two phases, this method runs the
   189     ///first phase. After the first phase the maximum flow value and a
   185     ///first phase. After the first phase the maximum flow value and a
   190     ///minimum value cut can already be computed, though a maximum flow
   186     ///minimum value cut can already be computed, though a maximum flow
   191     ///is net yet obtained. So after calling this method \ref flowValue
   187     ///is not yet obtained. So after calling this method \ref flowValue
   192     ///and \ref actMinCut gives proper results.
   188     ///and \ref actMinCut gives proper results.
   193     ///\warning: \ref minCut, \ref minMinCut and \ref maxMinCut do not
   189     ///\warning: \ref minCut, \ref minMinCut and \ref maxMinCut do not
   194     ///give minimum value cuts unless calling \ref preflowPhase2.
   190     ///give minimum value cuts unless calling \ref preflowPhase2.
   195     ///\pre The starting flow must be
   191     ///\pre The starting flow must be
   196     /// - a constant zero flow if \c fe is \c ZERO_FLOW,
   192     /// - a constant zero flow if \c fe is \c ZERO_FLOW,
   221       NNMap left(*g, INVALID);
   217       NNMap left(*g, INVALID);
   222       NNMap right(*g, INVALID);
   218       NNMap right(*g, INVALID);
   223       VecNode level_list(n,INVALID);
   219       VecNode level_list(n,INVALID);
   224       //List of the nodes in level i<n, set to n.
   220       //List of the nodes in level i<n, set to n.
   225 
   221 
   226       NodeIt v;
       
   227       for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
       
   228       //setting each node to level n
       
   229 
       
   230       if ( fe == NO_FLOW ) {
       
   231 	EdgeIt e;
       
   232 	for(g->first(e); g->valid(e); g->next(e)) flow->set(e,0);
       
   233       }
       
   234 
       
   235       switch (fe) { //computing the excess
       
   236       case PRE_FLOW:
       
   237 	{
       
   238 	  NodeIt v;
       
   239 	  for(g->first(v); g->valid(v); g->next(v)) {
       
   240 	    Num exc=0;
       
   241 
       
   242 	    InEdgeIt e;
       
   243 	    for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
       
   244 	    OutEdgeIt f;
       
   245 	    for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
       
   246 
       
   247 	    excess.set(v,exc);
       
   248 
       
   249 	    //putting the active nodes into the stack
       
   250 	    int lev=level[v];
       
   251 	    if ( exc > 0 && lev < n && v != t ) 
       
   252 	      {
       
   253 		next.set(v,first[lev]);
       
   254 		first[lev]=v;
       
   255 	      }
       
   256 	  }
       
   257 	  break;
       
   258 	}
       
   259       case GEN_FLOW:
       
   260 	{
       
   261 	  NodeIt v;
       
   262 	  for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
       
   263 
       
   264 	  Num exc=0;
       
   265 	  InEdgeIt e;
       
   266 	  for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
       
   267 	  OutEdgeIt f;
       
   268 	  for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
       
   269 	  excess.set(t,exc);
       
   270 	  break;
       
   271 	}
       
   272       case ZERO_FLOW:
       
   273       case NO_FLOW:
       
   274 	{
       
   275 	  NodeIt v;
       
   276 	  for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
       
   277 	  break;
       
   278 	}
       
   279       }
       
   280 
       
   281       preflowPreproc(fe, next, first, level_list, left, right);
   222       preflowPreproc(fe, next, first, level_list, left, right);
   282       //End of preprocessing
   223       //End of preprocessing
   283 
       
   284 
   224 
   285       //Push/relabel on the highest level active nodes.
   225       //Push/relabel on the highest level active nodes.
   286       while ( true ) {
   226       while ( true ) {
   287 	if ( b == 0 ) {
   227 	if ( b == 0 ) {
   288 	  if ( !what_heur && !end && k > 0 ) {
   228 	  if ( !what_heur && !end && k > 0 ) {
   392 	    level.set(w,++newlevel);
   332 	    level.set(w,++newlevel);
   393 	    next.set(w,first[newlevel]);
   333 	    next.set(w,first[newlevel]);
   394 	    first[newlevel]=w;
   334 	    first[newlevel]=w;
   395 	    b=newlevel;
   335 	    b=newlevel;
   396 	  }
   336 	  }
   397 	}  // if stack[b] is nonempty
   337 	} 
   398       } // while(true)
   338       } // while(true)
   399 
   339 
   400       status=AFTER_PRE_FLOW_PHASE_2;
   340       status=AFTER_PRE_FLOW_PHASE_2;
   401     }
   341     }
   402 
   342 
   411       for(InEdgeIt e(*g,t);g->valid(e);g->next(e)) a+=(*flow)[e];
   351       for(InEdgeIt e(*g,t);g->valid(e);g->next(e)) a+=(*flow)[e];
   412       for(OutEdgeIt e(*g,t);g->valid(e);g->next(e)) a-=(*flow)[e];
   352       for(OutEdgeIt e(*g,t);g->valid(e);g->next(e)) a-=(*flow)[e];
   413       return a;
   353       return a;
   414       //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan   
   354       //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan   
   415     }
   355     }
   416     Num flowValue2() const {
   356 
   417       return excess[t];
       
   418 //       Num a=0;
       
   419 //       for(InEdgeIt e(*g,t);g->valid(e);g->next(e)) a+=(*flow)[e];
       
   420 //       for(OutEdgeIt e(*g,t);g->valid(e);g->next(e)) a-=(*flow)[e];
       
   421 //       return a;
       
   422 //       //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan  
       
   423       
       
   424     }
       
   425 
   357 
   426     ///Returns a minimum value cut after calling \ref preflowPhase1.
   358     ///Returns a minimum value cut after calling \ref preflowPhase1.
   427 
   359 
   428     ///After the first phase of the preflow algorithm the maximum flow
   360     ///After the first phase of the preflow algorithm the maximum flow
   429     ///value and a minimum value cut can already be computed. This
   361     ///value and a minimum value cut can already be computed. This
   650 
   582 
   651       return newlevel;
   583       return newlevel;
   652     }
   584     }
   653 
   585 
   654 
   586 
       
   587 
   655     void preflowPreproc(FlowEnum fe, NNMap& next, VecFirst& first,
   588     void preflowPreproc(FlowEnum fe, NNMap& next, VecFirst& first,
   656 			VecNode& level_list, NNMap& left, NNMap& right)
   589 			VecNode& level_list, NNMap& left, NNMap& right)
   657     {
   590     {
       
   591       switch (fe) { //setting excess
       
   592 	case NO_FLOW: 
       
   593 	{
       
   594 	  EdgeIt e;
       
   595 	  for(g->first(e); g->valid(e); g->next(e)) flow->set(e,0);
       
   596 	  
       
   597 	  NodeIt v;
       
   598 	  for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
       
   599 	  break;
       
   600 	}
       
   601 	case ZERO_FLOW: 
       
   602 	{
       
   603 	  NodeIt v;
       
   604 	  for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
       
   605 	  break;
       
   606 	}
       
   607 	case GEN_FLOW:
       
   608 	{
       
   609 	  NodeIt v;
       
   610 	  for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
       
   611 
       
   612 	  Num exc=0;
       
   613 	  InEdgeIt e;
       
   614 	  for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
       
   615 	  OutEdgeIt f;
       
   616 	  for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
       
   617 	  excess.set(t,exc);
       
   618 	  break;
       
   619 	}
       
   620 	default: break;
       
   621       }
       
   622 
       
   623       NodeIt v;
       
   624       for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
       
   625       //setting each node to level n
       
   626       
   658       std::queue<Node> bfs_queue;
   627       std::queue<Node> bfs_queue;
   659 
   628 
       
   629 
   660       switch (fe) {
   630       switch (fe) {
   661       case NO_FLOW:   //flow is already set to const zero in this case
   631       case NO_FLOW:   //flow is already set to const zero
   662       case ZERO_FLOW:
   632       case ZERO_FLOW:
   663 	{
   633 	{
   664 	  //Reverse_bfs from t, to find the starting level.
   634 	  //Reverse_bfs from t, to find the starting level.
   665 	  level.set(t,0);
   635 	  level.set(t,0);
   666 	  bfs_queue.push(t);
   636 	  bfs_queue.push(t);
   691 	    {
   661 	    {
   692 	      Num c=(*capacity)[e];
   662 	      Num c=(*capacity)[e];
   693 	      if ( c <= 0 ) continue;
   663 	      if ( c <= 0 ) continue;
   694 	      Node w=g->head(e);
   664 	      Node w=g->head(e);
   695 	      if ( level[w] < n ) {
   665 	      if ( level[w] < n ) {
   696 		if ( excess[w] <= 0 && w!=t ) 
   666 		if ( excess[w] <= 0 && w!=t ) //putting into the stack
   697 		  {
   667 		  { 
   698 		    next.set(w,first[level[w]]);
   668 		    next.set(w,first[level[w]]);
   699 		    first[level[w]]=w;
   669 		    first[level[w]]=w;
   700 		  }
   670 		  }
   701 		flow->set(e, c);
   671 		flow->set(e, c);
   702 		excess.set(w, excess[w]+c);
   672 		excess.set(w, excess[w]+c);
   704 	    }
   674 	    }
   705 	  break;
   675 	  break;
   706 	}
   676 	}
   707 
   677 
   708       case GEN_FLOW:
   678       case GEN_FLOW:
   709       case PRE_FLOW:
       
   710 	{
   679 	{
   711 	  //Reverse_bfs from t in the residual graph,
   680 	  //Reverse_bfs from t in the residual graph,
   712 	  //to find the starting level.
   681 	  //to find the starting level.
   713 	  level.set(t,0);
   682 	  level.set(t,0);
   714 	  bfs_queue.push(t);
   683 	  bfs_queue.push(t);
   746 		level.set(w, l);
   715 		level.set(w, l);
   747 	      }
   716 	      }
   748 	    }
   717 	    }
   749 	  }
   718 	  }
   750 
   719 
   751 
       
   752 	  //the starting flow
   720 	  //the starting flow
   753 	  OutEdgeIt e;
   721 	  OutEdgeIt e;
   754 	  for(g->first(e,s); g->valid(e); g->next(e))
   722 	  for(g->first(e,s); g->valid(e); g->next(e))
   755 	    {
   723 	    {
   756 	      Num rem=(*capacity)[e]-(*flow)[e];
   724 	      Num rem=(*capacity)[e]-(*flow)[e];
   757 	      if ( rem <= 0 ) continue;
   725 	      if ( rem <= 0 ) continue;
   758 	      Node w=g->head(e);
   726 	      Node w=g->head(e);
   759 	      if ( level[w] < n ) {
   727 	      if ( level[w] < n ) {
   760 		if ( excess[w] <= 0 && w!=t )
   728 		if ( excess[w] <= 0 && w!=t ) //putting into the stack
   761 		  {
   729 		  {
   762 		    next.set(w,first[level[w]]);
   730 		    next.set(w,first[level[w]]);
   763 		    first[level[w]]=w;
   731 		    first[level[w]]=w;
   764 		  }   
   732 		  }   
   765 		flow->set(e, (*capacity)[e]);
   733 		flow->set(e, (*capacity)[e]);
   775 	      if ( level[w] < n ) {
   743 	      if ( level[w] < n ) {
   776 		if ( excess[w] <= 0 && w!=t )
   744 		if ( excess[w] <= 0 && w!=t )
   777 		  {
   745 		  {
   778 		    next.set(w,first[level[w]]);
   746 		    next.set(w,first[level[w]]);
   779 		    first[level[w]]=w;
   747 		    first[level[w]]=w;
   780 		  }   
   748 		  }  
   781 		excess.set(w, excess[w]+(*flow)[f]);
   749 		excess.set(w, excess[w]+(*flow)[f]);
   782 		flow->set(f, 0);
   750 		flow->set(f, 0);
   783 	      }
   751 	      }
   784 	    }
   752 	    }
   785 	  break;
   753 	  break;
       
   754 	} //case GEN_FLOW
       
   755     
       
   756       case PRE_FLOW:
       
   757 	{
       
   758 	  //Reverse_bfs from t in the residual graph,
       
   759 	  //to find the starting level.
       
   760 	  level.set(t,0);
       
   761 	  bfs_queue.push(t);
       
   762 
       
   763 	  while (!bfs_queue.empty()) {
       
   764 
       
   765 	    Node v=bfs_queue.front();
       
   766 	    bfs_queue.pop();
       
   767 	    int l=level[v]+1;
       
   768 
       
   769 	    InEdgeIt e;
       
   770 	    for(g->first(e,v); g->valid(e); g->next(e)) {
       
   771 	      if ( (*capacity)[e] <= (*flow)[e] ) continue;
       
   772 	      Node w=g->tail(e);
       
   773 	      if ( level[w] == n && w != s ) {
       
   774 		bfs_queue.push(w);
       
   775 		Node z=level_list[l];
       
   776 		if ( g->valid(z) ) left.set(z,w);
       
   777 		right.set(w,z);
       
   778 		level_list[l]=w;
       
   779 		level.set(w, l);
       
   780 	      }
       
   781 	    }
       
   782 
       
   783 	    OutEdgeIt f;
       
   784 	    for(g->first(f,v); g->valid(f); g->next(f)) {
       
   785 	      if ( 0 >= (*flow)[f] ) continue;
       
   786 	      Node w=g->head(f);
       
   787 	      if ( level[w] == n && w != s ) {
       
   788 		bfs_queue.push(w);
       
   789 		Node z=level_list[l];
       
   790 		if ( g->valid(z) ) left.set(z,w);
       
   791 		right.set(w,z);
       
   792 		level_list[l]=w;
       
   793 		level.set(w, l);
       
   794 	      }
       
   795 	    }
       
   796 	  }
       
   797 
       
   798 
       
   799 	  //the starting flow
       
   800 	  OutEdgeIt e;
       
   801 	  for(g->first(e,s); g->valid(e); g->next(e))
       
   802 	    {
       
   803 	      Num rem=(*capacity)[e]-(*flow)[e];
       
   804 	      if ( rem <= 0 ) continue;
       
   805 	      Node w=g->head(e);
       
   806 	      if ( level[w] < n ) {
       
   807 		flow->set(e, (*capacity)[e]);
       
   808 		excess.set(w, excess[w]+rem);
       
   809 	      }
       
   810 	    }
       
   811 
       
   812 	  InEdgeIt f;
       
   813 	  for(g->first(f,s); g->valid(f); g->next(f))
       
   814 	    {
       
   815 	      if ( (*flow)[f] <= 0 ) continue;
       
   816 	      Node w=g->tail(f);
       
   817 	      if ( level[w] < n ) {
       
   818 		excess.set(w, excess[w]+(*flow)[f]);
       
   819 		flow->set(f, 0);
       
   820 	      }
       
   821 	    }
       
   822 	  
       
   823 	  NodeIt w; //computing the excess
       
   824 	  for(g->first(w); g->valid(w); g->next(w)) {
       
   825 	    Num exc=0;
       
   826 
       
   827 	    InEdgeIt e;
       
   828 	    for(g->first(e,w); g->valid(e); g->next(e)) exc+=(*flow)[e];
       
   829 	    OutEdgeIt f;
       
   830 	    for(g->first(f,w); g->valid(f); g->next(f)) exc-=(*flow)[f];
       
   831 
       
   832 	    excess.set(w,exc);
       
   833 
       
   834 	    //putting the active nodes into the stack
       
   835 	    int lev=level[w];
       
   836 	    if ( exc > 0 && lev < n && w != t ) 
       
   837 	      {
       
   838 		next.set(w,first[lev]);
       
   839 		first[lev]=w;
       
   840 	      }
       
   841 	  }
       
   842 	  break;
   786 	} //case PRE_FLOW
   843 	} //case PRE_FLOW
   787       }
   844       }
   788     } //preflowPreproc
   845     } //preflowPreproc
   789 
       
   790 
   846 
   791 
   847 
   792     void relabel(Node w, int newlevel, NNMap& next, VecFirst& first,
   848     void relabel(Node w, int newlevel, NNMap& next, VecFirst& first,
   793 		 VecNode& level_list, NNMap& left,
   849 		 VecNode& level_list, NNMap& left,
   794 		 NNMap& right, int& b, int& k, bool what_heur )
   850 		 NNMap& right, int& b, int& k, bool what_heur )
   850 	  left.set(w,INVALID);
   906 	  left.set(w,INVALID);
   851 	  level_list[newlevel]=w;
   907 	  level_list[newlevel]=w;
   852 	}
   908 	}
   853       }
   909       }
   854     } //relabel
   910     } //relabel
       
   911 
       
   912     void printexcess() {////
       
   913       std::cout << "Excesses:" <<std::endl;
       
   914 
       
   915       NodeIt v;
       
   916       for(g->first(v); g->valid(v); g->next(v)) {
       
   917 	std::cout << 1+(g->id(v)) << ":" << excess[v]<<std::endl; 
       
   918       }
       
   919     }
       
   920 
       
   921  void printlevel() {////
       
   922       std::cout << "Levels:" <<std::endl;
       
   923 
       
   924       NodeIt v;
       
   925       for(g->first(v); g->valid(v); g->next(v)) {
       
   926 	std::cout << 1+(g->id(v)) << ":" << level[v]<<std::endl; 
       
   927       }
       
   928     }
       
   929 
       
   930 void printactive() {////
       
   931       std::cout << "Levels:" <<std::endl;
       
   932 
       
   933       NodeIt v;
       
   934       for(g->first(v); g->valid(v); g->next(v)) {
       
   935 	std::cout << 1+(g->id(v)) << ":" << level[v]<<std::endl; 
       
   936       }
       
   937     }
       
   938 
       
   939 
   855   };  //class MaxFlow
   940   };  //class MaxFlow
   856 } //namespace hugo
   941 } //namespace hugo
   857 
   942 
   858 #endif //HUGO_MAX_FLOW_H
   943 #endif //HUGO_MAX_FLOW_H
   859 
   944