[Lemon-commits] [lemon_svn] athos: r862 - hugo/trunk/src/work/athos

Lemon SVN svn at lemon.cs.elte.hu
Mon Nov 6 20:41:55 CET 2006


Author: athos
Date: Tue May 25 14:31:18 2004
New Revision: 862

Added:
   hugo/trunk/src/work/athos/min_cost_flow.cc
      - copied, changed from r859, /hugo/trunk/src/test/mincostflows_test.cc
Modified:
   hugo/trunk/src/work/athos/mincostflow.h

Log:
Almost ready.

Copied: hugo/trunk/src/work/athos/min_cost_flow.cc (from r859, /hugo/trunk/src/test/mincostflows_test.cc)
==============================================================================
--- /hugo/trunk/src/test/mincostflows_test.cc	(original)
+++ hugo/trunk/src/work/athos/min_cost_flow.cc	Tue May 25 14:31:18 2004
@@ -1,7 +1,7 @@
 #include <iostream>
 #include "test_tools.h"
 #include <hugo/list_graph.h>
-#include <hugo/mincostflows.h>
+#include <mincostflow.h>
 //#include <path.h>
 //#include <maps.h>
 
@@ -62,6 +62,7 @@
   length.set(v4_t, 8);
   length.set(v5_t, 8);
 
+  /*
   ListGraph::EdgeMap<int> capacity(graph);
 
   capacity.set(s_v1, 2);
@@ -72,36 +73,38 @@
   capacity.set(v3_v5, 1);
   capacity.set(v4_t, 1);
   capacity.set(v5_t, 2);
+  */
 
   //  ConstMap<Edge, int> const1map(1);
-  std::cout << "Mincostflows algorithm test..." << std::endl;
+  std::cout << "Enhanced capacity scaling algorithm test (for the mincostflow problem)..." << std::endl;
 
-  MinCostFlows< ListGraph, ListGraph::EdgeMap<int>, ListGraph::EdgeMap<int> >
-    surb_test(graph, length, capacity);
+  MinCostFlow< ListGraph, ListGraph::EdgeMap<int>, ListGraph::NodeMap<int> >
+    min_cost_flow_test(graph, length, supply_demand);
 
   int k=1;
 
-  check(  surb_test.run(s,t,k) == 1 && surb_test.totalLength() == 19,"One path, total length should be 19");
+  /*
+  check(  min_cost_flow_test.run(s,t,k) == 1 && min_cost_flow_test.totalLength() == 19,"One path, total length should be 19");
 
-  check(surb_test.checkComplementarySlackness(), "Is the primal-dual solution pair really optimal?");
+  check(min_cost_flow_test.checkComplementarySlackness(), "Is the primal-dual solution pair really optimal?");
   
   k=2;
   
-  check(  surb_test.run(s,t,k) == 2 && surb_test.totalLength() == 41,"Two paths, total length should be 41");
+  check(  min_cost_flow_test.run(s,t,k) == 2 && min_cost_flow_test.totalLength() == 41,"Two paths, total length should be 41");
 
-  check(surb_test.checkComplementarySlackness(), "Is the primal-dual solution pair really optimal?");
+  check(min_cost_flow_test.checkComplementarySlackness(), "Is the primal-dual solution pair really optimal?");
   
   
   k=4;
 
-  check(  surb_test.run(s,t,k) == 3 && surb_test.totalLength() == 64,"Three paths, total length should be 64");
+  check(  min_cost_flow_test.run(s,t,k) == 3 && min_cost_flow_test.totalLength() == 64,"Three paths, total length should be 64");
 
-  check(surb_test.checkComplementarySlackness(), "Is the primal-dual solution pair really optimal?");
+  check(min_cost_flow_test.checkComplementarySlackness(), "Is the primal-dual solution pair really optimal?");
 
 
   cout << (passed ? "All tests passed." : "Some of the tests failed!!!")
        << endl;
 
   return passed ? 0 : 1;
-
+  */
 }

Modified: hugo/trunk/src/work/athos/mincostflow.h
==============================================================================
--- hugo/trunk/src/work/athos/mincostflow.h	(original)
+++ hugo/trunk/src/work/athos/mincostflow.h	Tue May 25 14:31:18 2004
@@ -59,7 +59,7 @@
     class ModLengthMap {   
       //typedef typename ResGraphType::template NodeMap<Length> NodeMap;
       typedef typename Graph::template NodeMap<Length> NodeMap;
-      const ResGraphType& G;
+      const ResGraphType& res_graph;
       //      const EdgeIntMap& rev;
       const LengthMap &ol;
       const NodeMap &pot;
@@ -68,22 +68,22 @@
       typedef typename LengthMap::ValueType ValueType;
 	
       ValueType operator[](typename ResGraphType::Edge e) const {     
-	if (G.forward(e))
-	  return  ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);   
+	if (res_graph.forward(e))
+	  return  ol[e]-(pot[res_graph.head(e)]-pot[res_graph.tail(e)]);   
 	else
-	  return -ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);   
+	  return -ol[e]-(pot[res_graph.head(e)]-pot[res_graph.tail(e)]);   
       }     
 	
-      ModLengthMap(const ResGraphType& _G,
+      ModLengthMap(const ResGraphType& _res_graph,
 		   const LengthMap &o,  const NodeMap &p) : 
-	G(_G), /*rev(_rev),*/ ol(o), pot(p){}; 
+	res_graph(_res_graph), /*rev(_rev),*/ ol(o), pot(p){}; 
     };//ModLengthMap
 
 
   protected:
     
     //Input
-    const Graph& G;
+    const Graph& graph;
     const LengthMap& length;
     const SupplyDemandMap& supply_demand;//supply or demand of nodes
 
@@ -104,8 +104,8 @@
   public :
 
 
-    MinCostFlow(Graph& _G, LengthMap& _length, SupplyDemandMap& _supply_demand) : G(_G), 
-      length(_length), supply_demand(_supply_demand), flow(_G), potential(_G){ }
+    MinCostFlow(Graph& _graph, LengthMap& _length, SupplyDemandMap& _supply_demand) : graph(_graph), 
+      length(_length), supply_demand(_supply_demand), flow(_graph), potential(_graph){ }
 
     
     ///Runs the algorithm.
@@ -114,7 +114,7 @@
 
     ///\todo May be it does make sense to be able to start with a nonzero 
     /// feasible primal-dual solution pair as well.
-    int run() {
+    void run() {
 
       //Resetting variables from previous runs
       //total_length = 0;
@@ -124,17 +124,18 @@
 	std::greater<SupplyDemand> > 	HeapType;
 
       //A heap for the excess nodes
-      HeapMap excess_nodes_map(G,-1);
+      HeapMap excess_nodes_map(graph,-1);
       HeapType excess_nodes(excess_nodes_map);
 
       //A heap for the deficit nodes
-      HeapMap deficit_nodes_map(G,-1);
+      HeapMap deficit_nodes_map(graph,-1);
       HeapType deficit_nodes(deficit_nodes_map);
 
       //A container to store nonabundant arcs
       list<Edge> nonabundant_arcs;
-      
-      FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
+
+	
+      FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
 	flow.set(e,0);
 	nonabundant_arcs.push_back(e);
       }
@@ -145,12 +146,12 @@
       typedef UnionFindEnum<Node, Graph::template NodeMap> UFE;
 
       //A union-find structure to store the abundant components
-      UFE::MapType abund_comp_map(G);
+      UFE::MapType abund_comp_map(graph);
       UFE abundant_components(abund_comp_map);
 
 
 
-      FOR_EACH_LOC(typename Graph::NodeIt, n, G){
+      FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
        	excess_deficit.set(n,supply_demand[n]);
 	//A supply node
 	if (excess_deficit[n] > 0){
@@ -174,10 +175,39 @@
       //can be the maximum deficit here
       SupplyDemand max_excess = delta;
       
-      //We need a residual graph which is uncapacitated
-      ResGraphType res_graph(G, flow);
+      //ConstMap<Edge,SupplyDemand> ConstEdgeMap;
 
+      //We need a residual graph which is uncapacitated
+      ResGraphType res_graph(graph, flow);
+      
+      //An EdgeMap to tell which arcs are abundant
+      template typename Graph::EdgeMap<bool> abundant_arcs(graph);
 
+      //Let's construct the sugraph consisting only of the abundant edges
+      typedef ConstMap< typename Graph::Node, bool > ConstNodeMap;
+      ConstNodeMap const_true_map(true);
+      typedef SubGraphWrapper< Graph, ConstNodeMap, 
+	 template typename Graph::EdgeMap<bool> > 
+	AbundantGraph;
+      AbundantGraph abundant_graph(graph, const_true_map, abundant_arcs );
+      
+      //Let's construct the residual graph for the abundant graph
+      typedef ResGraphWrapper<const AbundantGraph,int,CapacityMap,EdgeIntMap> 
+	ResAbGraph;
+      //Again uncapacitated
+      ResAbGraph res_ab_graph(abundant_graph, flow);
+      
+      //We need things for the bfs
+      typename ResAbGraph::NodeMap<bool> bfs_reached(res_ab_graph);
+      typename ResAbGraph::NodeMap<typename ResAbGraph::Edge> 
+	bfs_pred(res_ab_graph); 
+      NullMap<typename ResAbGraph::Node, int> bfs_dist_dummy(res_ab_graph);
+      //We want to run bfs-es (more) on this graph 'res_ab_graph'
+      Bfs < ResAbGraph , 
+	typename ResAbGraph::NodeMap<bool>, 
+	typename ResAbGraph::NodeMap<typename ResAbGraph::Edge>,
+	NullMap<typename ResAbGraph::Node, int> > 
+	bfs(res_ab_graph, bfs_reached, bfs_pred, bfs_dist_dummy);
       
       ModLengthMap mod_length(res_graph, length, potential);
 
@@ -205,13 +235,55 @@
 	      Node b = abundant_components.find(res_graph.tail(i));
 	      //Merge
 	      if (a != b){
-		//Find path and augment
-		//!!!
-		//Find path and augment
 		abundant_components.join(a,b);
+		//We want to push the smaller
+		//Which has greater absolut value excess/deficit
+		Node root=(abs(excess_deficit[a])>abs(excess_deficit[b]))?a:b;
+		//Which is the other
+		Node non_root = ( a == root ) ? b : a ;
+		abundant_components.makeRep(root);
+		SupplyDemand qty_to_augment = abs(excess_deficit[non_root]); 
+		//Push the positive value
+		if (excess_deficit[non_root] < 0)
+		  swap(root, non_root);
+		//If the non_root node has excess/deficit at all
+		if (qty_to_augment>0){
+		  //Find path and augment
+		  bfs.run(non_root);
+		  //root should be reached
+		  
+		  //Augmenting on the found path
+		  Node n=root;
+		  ResGraphEdge e;
+		  while (n!=non_root){
+		    e = bfs_pred(n);
+		    n = res_graph.tail(e);
+		    res_graph.augment(e,qty_to_augment);
+		  }
+	  
+		  //We know that non_root had positive excess
+		  excess_nodes[non_root] -= qty_to_augment;
+		  //But what about root node
+		  //It might have been positive and so became larger
+		  if (excess_deficit[root]>0){
+		    excess_nodes[root] += qty_to_augment;
+		  }
+		  else{
+		    //Or negative but not turned into positive
+		    deficit_nodes[root] -= qty_to_augment;
+		  }
+
+		  //Update the excess_deficit map
+		  excess_deficit[non_root] -= qty_to_augment;
+		  excess_deficit[root] += qty_to_augment;
+
+		  
+		}
 	      }
 	      //What happens to i?
-	      nonabundant_arcs.erase(i);
+	      //Marci and Zsolt says I shouldn't do such things
+	      nonabundant_arcs.erase(i++);
+	      abundant_arcs[i] = true;
 	    }
 	    else
 	      ++i;
@@ -222,19 +294,19 @@
 	Node s = excess_nodes.top(); 
 	SupplyDemand max_excess = excess_nodes[s];
 	Node t = deficit_nodes.top(); 
-	if (max_excess < dificit_nodes[t]){
-	  max_excess = dificit_nodes[t];
+	if (max_excess < deficit_nodes[t]){
+	  max_excess = deficit_nodes[t];
 	}
 
 
-	while(max_excess > 0){
-
+	while(max_excess > (n-1)*delta/n){
+	  
 	  
 	  //s es t valasztasa
-
+	  
 	  //Dijkstra part	
 	  dijkstra.run(s);
-
+	  
 	  /*We know from theory that t can be reached
 	  if (!dijkstra.reached(t)){
 	    //There are no k paths from s to t
@@ -263,6 +335,11 @@
 	      total_length -= length[e];	    
 	    */
 	  }
+	  
+	  //Update the excess_deficit map
+	  excess_deficit[s] -= delta;
+	  excess_deficit[t] += delta;
+	  
 
 	  //Update the excess_nodes heap
 	  if (delta >= excess_nodes[s]){
@@ -285,6 +362,15 @@
 	    deficit_nodes[t] -= delta;
 	  }
 	  //Dijkstra part ends here
+	  
+	  //Choose s and t again
+	  s = excess_nodes.top(); 
+	  max_excess = excess_nodes[s];
+	  t = deficit_nodes.top(); 
+	  if (max_excess < deficit_nodes[t]){
+	    max_excess = deficit_nodes[t];
+	  }
+
 	}
 
 	/*
@@ -297,7 +383,7 @@
 	/*This is not necessary here
 	//Update the max_excess
 	max_excess = 0;
-	FOR_EACH_LOC(typename Graph::NodeIt, n, G){
+	FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
 	  if (max_excess < excess_deficit[n]){
 	    max_excess = excess_deficit[n];
 	  }
@@ -336,9 +422,9 @@
     bool checkComplementarySlackness(){
       Length mod_pot;
       Length fl_e;
-      FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
+      FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
 	//C^{\Pi}_{i,j}
-	mod_pot = length[e]-potential[G.head(e)]+potential[G.tail(e)];
+	mod_pot = length[e]-potential[graph.head(e)]+potential[graph.tail(e)];
 	fl_e = flow[e];
 	//	std::cout << fl_e << std::endl;
 	if (0<fl_e && fl_e<capacity[e]){



More information about the Lemon-commits mailing list