src/work/athos/mincostflow.h
changeset 658 b3564d0e9c60
parent 645 d93d8b9906d1
child 659 c5984e925384
equal deleted inserted replaced
2:0c59a6373a0b 3:97c30c04b2d4
     8 
     8 
     9 #include <hugo/dijkstra.h>
     9 #include <hugo/dijkstra.h>
    10 #include <hugo/graph_wrapper.h>
    10 #include <hugo/graph_wrapper.h>
    11 #include <hugo/maps.h>
    11 #include <hugo/maps.h>
    12 #include <vector>
    12 #include <vector>
       
    13 #include <list>
    13 #include <for_each_macros.h>
    14 #include <for_each_macros.h>
       
    15 #include <hugo/union_find.h>
    14 
    16 
    15 namespace hugo {
    17 namespace hugo {
    16 
    18 
    17 /// \addtogroup galgs
    19 /// \addtogroup galgs
    18 /// @{
    20 /// @{
   116 
   118 
   117       //Resetting variables from previous runs
   119       //Resetting variables from previous runs
   118       //total_length = 0;
   120       //total_length = 0;
   119 
   121 
   120       typedef typename Graph::template NodeMap<int> HeapMap;
   122       typedef typename Graph::template NodeMap<int> HeapMap;
   121       typedef Heap<Node, SupplyDemand, typename Graph::template NodeMap<int>,
   123       typedef Heap< Node, SupplyDemand, typename Graph::template NodeMap<int>,
   122 	std::greater<SupplyDemand> > 	HeapType;
   124 	std::greater<SupplyDemand> > 	HeapType;
   123 
   125 
   124       //A heap for the excess nodes
   126       //A heap for the excess nodes
   125       HeapMap excess_nodes_map(G,-1);
   127       HeapMap excess_nodes_map(G,-1);
   126       HeapType excess_nodes(excess_nodes_map);
   128       HeapType excess_nodes(excess_nodes_map);
   127 
   129 
   128       //A heap for the deficit nodes
   130       //A heap for the deficit nodes
   129       HeapMap deficit_nodes_map(G,-1);
   131       HeapMap deficit_nodes_map(G,-1);
   130       HeapType deficit_nodes(deficit_nodes_map);
   132       HeapType deficit_nodes(deficit_nodes_map);
   131 
   133 
       
   134       //A container to store nonabundant arcs
       
   135       list<Edge> nonabundant_arcs;
   132       
   136       
   133       FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
   137       FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
   134 	flow.set(e,0);
   138 	flow.set(e,0);
       
   139 	nonabundant_arcs.push_back(e);
   135       }
   140       }
   136 
   141 
   137       //Initial value for delta
   142       //Initial value for delta
   138       SupplyDemand delta = 0;
   143       SupplyDemand delta = 0;
       
   144 
       
   145       typedef UnionFindEnum<Node, Graph::template NodeMap> UFE;
       
   146 
       
   147       //A union-find structure to store the abundant components
       
   148       UFE::MapType abund_comp_map(G);
       
   149       UFE abundant_components(abund_comp_map);
       
   150 
       
   151 
   139 
   152 
   140       FOR_EACH_LOC(typename Graph::NodeIt, n, G){
   153       FOR_EACH_LOC(typename Graph::NodeIt, n, G){
   141        	excess_deficit.set(n,supply_demand[n]);
   154        	excess_deficit.set(n,supply_demand[n]);
   142 	//A supply node
   155 	//A supply node
   143 	if (excess_deficit[n] > 0){
   156 	if (excess_deficit[n] > 0){
   151 	if (delta < abs(excess_deficit[n])){
   164 	if (delta < abs(excess_deficit[n])){
   152 	  delta = abs(excess_deficit[n]);
   165 	  delta = abs(excess_deficit[n]);
   153 	}
   166 	}
   154 	//Initialize the copy of the Dijkstra potential to zero
   167 	//Initialize the copy of the Dijkstra potential to zero
   155 	potential.set(n,0);
   168 	potential.set(n,0);
       
   169 	//Every single point is an abundant component initially 
       
   170 	abundant_components.insert(n);
   156       }
   171       }
   157 
   172 
   158       //It'll be allright as an initial value, though this value 
   173       //It'll be allright as an initial value, though this value 
   159       //can be the maximum deficit here
   174       //can be the maximum deficit here
   160       SupplyDemand max_excess = delta;
   175       SupplyDemand max_excess = delta;
   169       Dijkstra<ResGraphType, ModLengthMap> dijkstra(res_graph, mod_length);
   184       Dijkstra<ResGraphType, ModLengthMap> dijkstra(res_graph, mod_length);
   170 
   185 
   171 
   186 
   172       while (max_excess > 0){
   187       while (max_excess > 0){
   173 
   188 
       
   189 	//Reset delta if still too big
       
   190 	if (8*number_of_nodes*max_excess <= delta){
       
   191 	  delta = max_excess;
       
   192 	  
       
   193 	}
       
   194 
   174 	/*
   195 	/*
   175 	 * Beginning of the delta scaling phase 
   196 	 * Beginning of the delta scaling phase 
   176 	*/
   197 	*/
   177 	
       
   178 	//Merge and stuff
   198 	//Merge and stuff
       
   199 	{
       
   200 	  SupplyDemand buf=8*number_of_nodes*delta;
       
   201 	  list<Edge>::iterator i = nonabundant_arcs.begin();
       
   202 	  while ( i != nonabundant_arcs.end() ){
       
   203 	    if (flow[i]>=buf){
       
   204 	      Node a = abundant_components.find(res_graph.head(i));
       
   205 	      Node b = abundant_components.find(res_graph.tail(i));
       
   206 	      //Merge
       
   207 	      if (a != b){
       
   208 		//Find path and augment
       
   209 		//!!!
       
   210 		//Find path and augment
       
   211 		abundant_components.join(a,b);
       
   212 	      }
       
   213 	      //What happens to i?
       
   214 	      nonabundant_arcs.erase(i);
       
   215 	    }
       
   216 	    else
       
   217 	      ++i;
       
   218 	  }
       
   219 	}
       
   220 
   179 
   221 
   180 	Node s = excess_nodes.top(); 
   222 	Node s = excess_nodes.top(); 
   181 	SupplyDemand max_excess = excess_nodes[s];
   223 	SupplyDemand max_excess = excess_nodes[s];
   182 	Node t = deficit_nodes.top(); 
   224 	Node t = deficit_nodes.top(); 
   183 	if (max_excess < dificit_nodes[t]){
   225 	if (max_excess < dificit_nodes[t]){
   259 	  if (max_excess < excess_deficit[n]){
   301 	  if (max_excess < excess_deficit[n]){
   260 	    max_excess = excess_deficit[n];
   302 	    max_excess = excess_deficit[n];
   261 	  }
   303 	  }
   262 	}
   304 	}
   263 	*/
   305 	*/
   264 	//Reset delta if still too big
   306 
   265 	if (8*number_of_nodes*max_excess <= delta){
       
   266 	  delta = max_excess;
       
   267 	  
       
   268 	}
       
   269 	  
   307 	  
   270       }//while(max_excess > 0)
   308       }//while(max_excess > 0)
   271       
   309       
   272 
   310 
   273       return i;
   311       return i;