# HG changeset patch
# User athos
# Date 1085488278 0
# Node ID c5984e925384d09a523cfbc5c4481188abdf11b2
# Parent  b3564d0e9c60ff31f1fd7b1a8f9bd2a68d770757
Almost ready.

diff -r b3564d0e9c60 -r c5984e925384 src/work/athos/min_cost_flow.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/athos/min_cost_flow.cc	Tue May 25 12:31:18 2004 +0000
@@ -0,0 +1,110 @@
+#include <iostream>
+#include "test_tools.h"
+#include <hugo/list_graph.h>
+#include <mincostflow.h>
+//#include <path.h>
+//#include <maps.h>
+
+using namespace std;
+using namespace hugo;
+
+
+
+bool passed = true;
+/*
+void check(bool rc, char *msg="") {
+  passed = passed && rc;
+  if(!rc) {
+    std::cerr << "Test failed! ("<< msg << ")" << std::endl; \
+ 
+
+  }
+}
+*/
+
+
+int main()
+{
+
+  typedef ListGraph::Node Node;
+  typedef ListGraph::Edge Edge;
+
+  ListGraph graph;
+
+  //Ahuja könyv példája
+
+  Node s=graph.addNode();
+  Node v1=graph.addNode();  
+  Node v2=graph.addNode();
+  Node v3=graph.addNode();
+  Node v4=graph.addNode();
+  Node v5=graph.addNode();
+  Node t=graph.addNode();
+
+  Edge s_v1=graph.addEdge(s, v1);
+  Edge v1_v2=graph.addEdge(v1, v2);
+  Edge s_v3=graph.addEdge(s, v3);
+  Edge v2_v4=graph.addEdge(v2, v4);
+  Edge v2_v5=graph.addEdge(v2, v5);
+  Edge v3_v5=graph.addEdge(v3, v5);
+  Edge v4_t=graph.addEdge(v4, t);
+  Edge v5_t=graph.addEdge(v5, t);
+  
+
+  ListGraph::EdgeMap<int> length(graph);
+
+  length.set(s_v1, 6);
+  length.set(v1_v2, 4);
+  length.set(s_v3, 10);
+  length.set(v2_v4, 5);
+  length.set(v2_v5, 1);
+  length.set(v3_v5, 4);
+  length.set(v4_t, 8);
+  length.set(v5_t, 8);
+
+  /*
+  ListGraph::EdgeMap<int> capacity(graph);
+
+  capacity.set(s_v1, 2);
+  capacity.set(v1_v2, 2);
+  capacity.set(s_v3, 1);
+  capacity.set(v2_v4, 1);
+  capacity.set(v2_v5, 1);
+  capacity.set(v3_v5, 1);
+  capacity.set(v4_t, 1);
+  capacity.set(v5_t, 2);
+  */
+
+  //  ConstMap<Edge, int> const1map(1);
+  std::cout << "Enhanced capacity scaling algorithm test (for the mincostflow problem)..." << std::endl;
+
+  MinCostFlow< ListGraph, ListGraph::EdgeMap<int>, ListGraph::NodeMap<int> >
+    min_cost_flow_test(graph, length, supply_demand);
+
+  int k=1;
+
+  /*
+  check(  min_cost_flow_test.run(s,t,k) == 1 && min_cost_flow_test.totalLength() == 19,"One path, total length should be 19");
+
+  check(min_cost_flow_test.checkComplementarySlackness(), "Is the primal-dual solution pair really optimal?");
+  
+  k=2;
+  
+  check(  min_cost_flow_test.run(s,t,k) == 2 && min_cost_flow_test.totalLength() == 41,"Two paths, total length should be 41");
+
+  check(min_cost_flow_test.checkComplementarySlackness(), "Is the primal-dual solution pair really optimal?");
+  
+  
+  k=4;
+
+  check(  min_cost_flow_test.run(s,t,k) == 3 && min_cost_flow_test.totalLength() == 64,"Three paths, total length should be 64");
+
+  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;
+  */
+}
diff -r b3564d0e9c60 -r c5984e925384 src/work/athos/mincostflow.h
--- a/src/work/athos/mincostflow.h	Mon May 24 14:13:03 2004 +0000
+++ b/src/work/athos/mincostflow.h	Tue May 25 12:31:18 2004 +0000
@@ -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;
       
+      //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);
 
@@ -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]){