Index: src/work/athos/mincostflow.h
===================================================================
--- src/work/athos/mincostflow.h	(revision 657)
+++ src/work/athos/mincostflow.h	(revision 659)
@@ -60,5 +60,5 @@
       //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;
@@ -69,13 +69,13 @@
 	
       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
 
@@ -84,5 +84,5 @@
     
     //Input
-    const Graph& G;
+    const Graph& graph;
     const LengthMap& length;
     const SupplyDemandMap& supply_demand;//supply or demand of nodes
@@ -105,6 +105,6 @@
 
 
-    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){ }
 
     
@@ -115,5 +115,5 @@
     ///\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
@@ -125,15 +125,16 @@
 
       //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);
@@ -146,10 +147,10 @@
 
       //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
@@ -175,8 +176,37 @@
       SupplyDemand max_excess = delta;
       
+      //ConstMap<Edge,SupplyDemand> ConstEdgeMap;
+
       //We need a residual graph which is uncapacitated
-      ResGraphType res_graph(G, flow);
-
-
+      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);
@@ -206,11 +236,53 @@
 	      //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
@@ -223,17 +295,17 @@
 	SupplyDemand max_excess = excess_nodes[s];
 	Node t = deficit_nodes.top(); 
-	if (max_excess < dificit_nodes[t]){
-	  max_excess = dificit_nodes[t];
-	}
-
-
-	while(max_excess > 0){
-
+	if (max_excess < deficit_nodes[t]){
+	  max_excess = deficit_nodes[t];
+	}
+
+
+	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)){
@@ -264,4 +336,9 @@
 	    */
 	  }
+	  
+	  //Update the excess_deficit map
+	  excess_deficit[s] -= delta;
+	  excess_deficit[t] += delta;
+	  
 
 	  //Update the excess_nodes heap
@@ -286,4 +363,13 @@
 	  }
 	  //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];
+	  }
+
 	}
 
@@ -298,5 +384,5 @@
 	//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];
@@ -337,7 +423,7 @@
       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;
