# HG changeset patch
# User athos
# Date 1084470143 0
# Node ID 933f593824c24740078a3aed562fb762995c916f
# Parent  aacabcd724f0d99de5fa3858495e8aa02d2ab384
Started mincostflow.

diff -r aacabcd724f0 -r 933f593824c2 src/work/athos/mincostflow.h
--- a/src/work/athos/mincostflow.h	Thu May 13 17:33:40 2004 +0000
+++ b/src/work/athos/mincostflow.h	Thu May 13 17:42:23 2004 +0000
@@ -36,13 +36,13 @@
   /// where \c M is the value of the maximal flow.
   ///
   ///\author Attila Bernath
-  template <typename Graph, typename LengthMap, typename SupplyMap>
+  template <typename Graph, typename LengthMap, typename SupplyDemandMap>
   class MinCostFlow {
 
     typedef typename LengthMap::ValueType Length;
 
 
-    typedef typename SupplyMap::ValueType Supply;
+    typedef typename SupplyDemandMap::ValueType SupplyDemand;
     
     typedef typename Graph::Node Node;
     typedef typename Graph::NodeIt NodeIt;
@@ -84,7 +84,7 @@
     //Input
     const Graph& G;
     const LengthMap& length;
-    const SupplyMap& supply;//supply or demand of nodes
+    const SupplyDemandMap& supply_demand;//supply or demand of nodes
 
 
     //auxiliary variables
@@ -94,7 +94,7 @@
     //To store the potentila (dual variables)
     typename Graph::template NodeMap<Length> potential;
     //To store excess-deficit values
-    SupplyMap excess;
+    SupplyDemandMap excess_deficit;
     
 
     Length total_length;
@@ -103,39 +103,62 @@
   public :
 
 
-    MinCostFlow(Graph& _G, LengthMap& _length, SupplyMap& _supply) : G(_G), 
-      length(_length), supply(_supply), flow(_G), potential(_G){ }
+    MinCostFlow(Graph& _G, LengthMap& _length, SupplyDemandMap& _supply_demand) : G(_G), 
+      length(_length), supply_demand(_supply_demand), flow(_G), potential(_G){ }
 
     
     ///Runs the algorithm.
 
     ///Runs the algorithm.
-    ///Returns k if there are at least k edge-disjoint paths from s to t.
-    ///Otherwise it returns the number of found edge-disjoint paths from s to t.
+
     ///\todo May be it does make sense to be able to start with a nonzero 
     /// feasible primal-dual solution pair as well.
     int run() {
 
       //Resetting variables from previous runs
-      total_length = 0;
+      //total_length = 0;
+
+      typedef typename Graph::template NodeMap<int> HeapMap;
+      typedef Heap<Node, SupplyDemand, typename Graph::template NodeMap<int>,
+	std::greater<SupplyDemand> > 	HeapType;
+
+      //A heap for the excess nodes
+      HeapMap excess_nodes_map(G,-1);
+      HeapType excess_nodes(excess_nodes_map);
+
+      //A heap for the deficit nodes
+      HeapMap deficit_nodes_map(G,-1);
+      HeapType deficit_nodes(deficit_nodes_map);
+
       
       FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
 	flow.set(e,0);
       }
 
       //Initial value for delta
-      Supply delta = 0;
-      
+      SupplyDemand delta = 0;
+
       FOR_EACH_LOC(typename Graph::NodeIt, n, G){
-	if (delta < abs(supply[e])){
-	  delta = abs(supply[e]);
+       	excess_deficit.set(n,supply_demand[n]);
+	//A supply node
+	if (excess_deficit[n] > 0){
+	  excess_nodes.push(n,excess_deficit[n]);
 	}
-	excess.set(n,supply[e]);
+	//A demand node
+	if (excess_deficit[n] < 0){
+	  deficit_nodes.push(n, - excess_deficit[n]);
+	}
+	//Finding out starting value of delta
+	if (delta < abs(excess_deficit[n])){
+	  delta = abs(excess_deficit[n]);
+	}
 	//Initialize the copy of the Dijkstra potential to zero
 	potential.set(n,0);
       }
-      
 
+      //It'll be allright as an initial value, though this value 
+      //can be the maximum deficit here
+      SupplyDemand max_excess = delta;
       
       //We need a residual graph which is uncapacitated
       ResGraphType res_graph(G, flow);
@@ -147,45 +170,102 @@
       Dijkstra<ResGraphType, ModLengthMap> dijkstra(res_graph, mod_length);
 
 
-      int i;
-      for (i=0; i<k; ++i){
-	dijkstra.run(s);
-	if (!dijkstra.reached(t)){
-	  //There are no k paths from s to t
-	  break;
-	};
+      while (max_excess > 0){
+
 	
-	//We have to copy the potential
-	FOR_EACH_LOC(typename ResGraphType::NodeIt, n, res_graph){
-	  potential[n] += dijkstra.distMap()[n];
+	//Merge and stuff
+
+	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];
+	}
+
+
+	while(max_excess > ){
+
+	  
+	  //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
+	    break;
+	  };
+	  */
+	  
+	  //We have to change the potential
+	  FOR_EACH_LOC(typename ResGraphType::NodeIt, n, res_graph){
+	    potential[n] += dijkstra.distMap()[n];
+	  }
+
+
+	  //Augmenting on the sortest path
+	  Node n=t;
+	  ResGraphEdge e;
+	  while (n!=s){
+	    e = dijkstra.pred(n);
+	    n = dijkstra.predNode(n);
+	    res_graph.augment(e,delta);
+	    /*
+	    //Let's update the total length
+	    if (res_graph.forward(e))
+	      total_length += length[e];
+	    else 
+	      total_length -= length[e];	    
+	    */
+	  }
+
+	  //Update the excess_nodes heap
+	  if (delta >= excess_nodes[s]){
+	    if (delta > excess_nodes[s])
+	      deficit_nodes.push(s,delta - excess_nodes[s]);
+	    excess_nodes.pop();
+	    
+	  } 
+	  else{
+	    excess_nodes[s] -= delta;
+	  }
+	  //Update the deficit_nodes heap
+	  if (delta >= deficit_nodes[t]){
+	    if (delta > deficit_nodes[t])
+	      excess_nodes.push(t,delta - deficit_nodes[t]);
+	    deficit_nodes.pop();
+	    
+	  } 
+	  else{
+	    deficit_nodes[t] -= delta;
+	  }
+	  //Dijkstra part ends here
 	}
 
 	/*
-	{
+	 * End of the delta scaling phase 
+	*/
 
-	  typename ResGraphType::NodeIt n;
-	  for ( res_graph.first(n) ; res_graph.valid(n) ; res_graph.next(n) ) {
-	      potential[n] += dijkstra.distMap()[n];
+	//Whatever this means
+	delta = delta / 2;
+
+	/*This is not necessary here
+	//Update the max_excess
+	max_excess = 0;
+	FOR_EACH_LOC(typename Graph::NodeIt, n, G){
+	  if (max_excess < excess_deficit[n]){
+	    max_excess = excess_deficit[n];
 	  }
 	}
 	*/
-
-	//Augmenting on the sortest path
-	Node n=t;
-	ResGraphEdge e;
-	while (n!=s){
-	  e = dijkstra.pred(n);
-	  n = dijkstra.predNode(n);
-	  res_graph.augment(e,delta);
-	  //Let's update the total length
-	  if (res_graph.forward(e))
-	    total_length += length[e];
-	  else 
-	    total_length -= length[e];	    
+	//Reset delta if still too big
+	if (8*number_of_nodes*max_excess <= delta){
+	  delta = max_excess;
+	  
 	}
-
 	  
-      }
+      }//while(max_excess > 0)
       
 
       return i;