Changeset 659:c5984e925384 in lemon-0.x for src/work

Timestamp:

05/25/04 14:31:18 (21 years ago)

Author:

athos

Branch:

default

Phase:

public

Convert:

svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@862

Message:

Almost ready.

Location:

src/work/athos

Files:

• min_cost_flow.cc (added)
• mincostflow.h (modified) (14 diffs)

src/work/athos/mincostflow.h

@@ -60 +60 @@
       //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 +69 @@
         
       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 +84 @@
     
     //Input
-    const Graph& G;
+    const Graph& graph;
     const LengthMap& length;
     const SupplyDemandMap& supply_demand;//supply or demand of nodes
@@ -105 +105 @@
 
 
-    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 +115 @@
     ///\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 +125 @@
 
       //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 +147 @@
 
       //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 +176 @@
       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 +236 @@
               //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 +295 @@
         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 +336 @@
             */
           }
+          
+          //Update the excess_deficit map
+          excess_deficit[s] -= delta;
+          excess_deficit[t] += delta;
+          
 
           //Update the excess_nodes heap
@@ -286 +363 @@
           }
           //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 +384 @@
         //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 +423 @@
       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;