Changeset 2625:c51b320bc51c in lemon-0.x for lemon/cost_scaling.h

Timestamp:

10/15/08 14:04:11 (16 years ago)

Author:

Peter Kovacs

Branch:

default

Phase:

public

Convert:

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

Message:

Major improvement in the cost scaling algorithm

• Add a new variant that use the partial augment-relabel method.
• Use this method instead of push-relabel by default.
• Use the "Early Termination" heuristic instead of "Price Refinement".

Using the new method and heuristic the algorithm proved to be
2-2.5 times faster on all input files.

File:

• lemon/cost_scaling.h (modified) (22 diffs)

lemon/cost_scaling.h

@@ -21 +21 @@
 
 /// \ingroup min_cost_flow
-///
 /// \file
 /// \brief Cost scaling algorithm for finding a minimum cost flow.
@@ -43 +42 @@
   ///
   /// \ref CostScaling implements the cost scaling algorithm performing
-  /// generalized push-relabel operations for finding a minimum cost
+  /// augment/push and relabel operations for finding a minimum cost
   /// flow.
   ///
@@ -117 +116 @@
 
       ///\e
-      typename Map::Value operator[](const ResEdge &e) const {
-        return ResGraph::forward(e) ?  _cost_map[e] : -_cost_map[e];
+      inline typename Map::Value operator[](const ResEdge &e) const {
+        return ResGraph::forward(e) ? _cost_map[e] : -_cost_map[e];
       }
 
@@ -143 +142 @@
 
       ///\e
-      LCost operator[](const Edge &e) const {
+      inline LCost operator[](const Edge &e) const {
         return _cost_map[e] + _pot_map[_gr.source(e)]
                             - _pot_map[_gr.target(e)];
@@ -149 +148 @@
 
     }; //class ReducedCostMap
-
-  private:
-
-    // Scaling factor
-    static const int ALPHA = 4;
-
-    // Paramters for heuristics
-    static const int BF_HEURISTIC_EPSILON_BOUND = 5000;
-    static const int BF_HEURISTIC_BOUND_FACTOR  = 3;
 
   private:
@@ -192 +182 @@
     // The epsilon parameter used for cost scaling
     LCost _epsilon;
+    // The scaling factor
+    int _alpha;
 
   public:
@@ -214 +206 @@
       _res_graph(NULL), _red_cost(NULL), _excess(graph, 0)
     {
-      // Removing non-zero lower bounds
+      // Remove non-zero lower bounds
       _capacity = subMap(capacity, lower);
       Supply sum = 0;
@@ -246 +238 @@
       _res_graph(NULL), _red_cost(NULL), _excess(graph, 0)
     {
-      // Checking the sum of supply values
+      // Check the sum of supply values
       Supply sum = 0;
       for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
@@ -275 +267 @@
       _res_graph(NULL), _red_cost(NULL), _excess(graph, 0)
     {
-      // Removing nonzero lower bounds
+      // Remove nonzero lower bounds
       _capacity = subMap(capacity, lower);
       for (NodeIt n(_graph); n != INVALID; ++n) {
@@ -360 +352 @@
     /// Run the algorithm.
     ///
+    /// \param partial_augment By default the algorithm performs
+    /// partial augment and relabel operations in the cost scaling
+    /// phases. Set this parameter to \c false for using local push and
+    /// relabel operations instead.
+    ///
     /// \return \c true if a feasible flow can be found.
-    bool run() {
-      return init() && start();
+    bool run(bool partial_augment = true) {
+      if (partial_augment) {
+        return init() && startPartialAugment();
+      } else {
+        return init() && startPushRelabel();
+      }
     }
 
@@ -434 +435 @@
     bool init() {
       if (!_valid_supply) return false;
-
-      // Initializing flow and potential maps
+      // The scaling factor
+      _alpha = 8;
+
+      // Initialize flow and potential maps
       if (!_flow) {
         _flow = new FlowMap(_graph);
@@ -448 +451 @@
       _res_graph = new ResGraph(_graph, _capacity, *_flow);
 
-      // Initializing the scaled cost map and the epsilon parameter
+      // Initialize the scaled cost map and the epsilon parameter
       Cost max_cost = 0;
       int node_num = countNodes(_graph);
       for (EdgeIt e(_graph); e != INVALID; ++e) {
-        _cost[e] = LCost(_orig_cost[e]) * node_num * ALPHA;
+        _cost[e] = LCost(_orig_cost[e]) * node_num * _alpha;
         if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e];
       }
       _epsilon = max_cost * node_num;
 
-      // Finding a feasible flow using Circulation
+      // Find a feasible flow using Circulation
       Circulation< Graph, ConstMap<Edge, Capacity>, CapacityEdgeMap,
                    SupplyMap >
@@ -465 +468 @@
     }
 
-
-    /// Execute the algorithm.
-    bool start() {
+    /// Execute the algorithm performing partial augmentation and
+    /// relabel operations.
+    bool startPartialAugment() {
+      // Paramters for heuristics
+      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
+      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
+      // Maximum augment path length
+      const int MAX_PATH_LENGTH = 4;
+
+      // Variables
+      typename Graph::template NodeMap<Edge> pred_edge(_graph);
+      typename Graph::template NodeMap<bool> forward(_graph);
+      typename Graph::template NodeMap<OutEdgeIt> next_out(_graph);
+      typename Graph::template NodeMap<InEdgeIt> next_in(_graph);
+      typename Graph::template NodeMap<bool> next_dir(_graph);
       std::deque<Node> active_nodes;
-      typename Graph::template NodeMap<bool> hyper(_graph, false);
+      std::vector<Node> path_nodes;
 
       int node_num = countNodes(_graph);
-      for ( ; _epsilon >= 1; _epsilon = _epsilon < ALPHA && _epsilon > 1 ?
-                                        1 : _epsilon / ALPHA )
+      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
+                                        1 : _epsilon / _alpha )
       {
-        // Performing price refinement heuristic using Bellman-Ford
-        // algorithm
+        // "Early Termination" heuristic: use Bellman-Ford algorithm
+        // to check if the current flow is optimal
         if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
           typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
-          ShiftCostMap shift_cost(_res_cost, _epsilon);
+          ShiftCostMap shift_cost(_res_cost, 1);
           BellmanFord<ResGraph, ShiftCostMap> bf(*_res_graph, shift_cost);
           bf.init(0);
@@ -486 +501 @@
           for (int i = 0; i < K && !done; ++i)
             done = bf.processNextWeakRound();
-          if (done) {
-            for (NodeIt n(_graph); n != INVALID; ++n)
-              (*_potential)[n] = bf.dist(n);
-            continue;
-          }
-        }
-
-        // Saturating edges not satisfying the optimality condition
+          if (done) break;
+        }
+
+        // Saturate edges not satisfying the optimality condition
         Capacity delta;
         for (EdgeIt e(_graph); e != INVALID; ++e) {
@@ -509 +520 @@
         }
 
-        // Finding active nodes (i.e. nodes with positive excess)
-        for (NodeIt n(_graph); n != INVALID; ++n)
+        // Find active nodes (i.e. nodes with positive excess)
+        for (NodeIt n(_graph); n != INVALID; ++n) {
           if (_excess[n] > 0) active_nodes.push_back(n);
-
-        // Performing push and relabel operations
+        }
+
+        // Initialize the next edge maps
+        for (NodeIt n(_graph); n != INVALID; ++n) {
+          next_out[n] = OutEdgeIt(_graph, n);
+          next_in[n] = InEdgeIt(_graph, n);
+          next_dir[n] = true;
+        }
+
+        // Perform partial augment and relabel operations
         while (active_nodes.size() > 0) {
+          // Select an active node (FIFO selection)
+          if (_excess[active_nodes[0]] <= 0) {
+            active_nodes.pop_front();
+            continue;
+          }
+          Node start = active_nodes[0];
+          path_nodes.clear();
+          path_nodes.push_back(start);
+
+          // Find an augmenting path from the start node
+          Node u, tip = start;
+          LCost min_red_cost;
+          while ( _excess[tip] >= 0 &&
+                  int(path_nodes.size()) <= MAX_PATH_LENGTH )
+          {
+            if (next_dir[tip]) {
+              for (OutEdgeIt e = next_out[tip]; e != INVALID; ++e) {
+                if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
+                  u = _graph.target(e);
+                  pred_edge[u] = e;
+                  forward[u] = true;
+                  next_out[tip] = e;
+                  tip = u;
+                  path_nodes.push_back(tip);
+                  goto next_step;
+                }
+              }
+              next_dir[tip] = false;
+            }
+            for (InEdgeIt e = next_in[tip]; e != INVALID; ++e) {
+              if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
+                u = _graph.source(e);
+                pred_edge[u] = e;
+                forward[u] = false;
+                next_in[tip] = e;
+                tip = u;
+                path_nodes.push_back(tip);
+                goto next_step;
+              }
+            }
+
+            // Relabel tip node
+            min_red_cost = std::numeric_limits<LCost>::max() / 2;
+            for (OutEdgeIt oe(_graph, tip); oe != INVALID; ++oe) {
+              if ( _capacity[oe] - (*_flow)[oe] > 0 &&
+                   (*_red_cost)[oe] < min_red_cost )
+                min_red_cost = (*_red_cost)[oe];
+            }
+            for (InEdgeIt ie(_graph, tip); ie != INVALID; ++ie) {
+              if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
+                min_red_cost = -(*_red_cost)[ie];
+            }
+            (*_potential)[tip] -= min_red_cost + _epsilon;
+
+            // Reset the next edge maps
+            next_out[tip] = OutEdgeIt(_graph, tip);
+            next_in[tip] = InEdgeIt(_graph, tip);
+            next_dir[tip] = true;
+
+            // Step back
+            if (tip != start) {
+              path_nodes.pop_back();
+              tip = path_nodes[path_nodes.size()-1];
+            }
+
+          next_step:
+            continue;
+          }
+
+          // Augment along the found path (as much flow as possible)
+          Capacity delta;
+          for (int i = 1; i < int(path_nodes.size()); ++i) {
+            u = path_nodes[i];
+            delta = forward[u] ?
+              _capacity[pred_edge[u]] - (*_flow)[pred_edge[u]] :
+              (*_flow)[pred_edge[u]];
+            delta = std::min(delta, _excess[path_nodes[i-1]]);
+            (*_flow)[pred_edge[u]] += forward[u] ? delta : -delta;
+            _excess[path_nodes[i-1]] -= delta;
+            _excess[u] += delta;
+            if (_excess[u] > 0 && _excess[u] <= delta) active_nodes.push_back(u);
+          }
+        }
+      }
+
+      // Compute node potentials for the original costs
+      ResidualCostMap<CostMap> res_cost(_orig_cost);
+      BellmanFord< ResGraph, ResidualCostMap<CostMap> >
+        bf(*_res_graph, res_cost);
+      bf.init(0); bf.start();
+      for (NodeIt n(_graph); n != INVALID; ++n)
+        (*_potential)[n] = bf.dist(n);
+
+      // Handle non-zero lower bounds
+      if (_lower) {
+        for (EdgeIt e(_graph); e != INVALID; ++e)
+          (*_flow)[e] += (*_lower)[e];
+      }
+      return true;
+    }
+
+    /// Execute the algorithm performing push and relabel operations.
+    bool startPushRelabel() {
+      // Paramters for heuristics
+      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
+      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
+
+      typename Graph::template NodeMap<bool> hyper(_graph, false);
+      typename Graph::template NodeMap<Edge> pred_edge(_graph);
+      typename Graph::template NodeMap<bool> forward(_graph);
+      typename Graph::template NodeMap<OutEdgeIt> next_out(_graph);
+      typename Graph::template NodeMap<InEdgeIt> next_in(_graph);
+      typename Graph::template NodeMap<bool> next_dir(_graph);
+      std::deque<Node> active_nodes;
+
+      int node_num = countNodes(_graph);
+      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
+                                        1 : _epsilon / _alpha )
+      {
+        // "Early Termination" heuristic: use Bellman-Ford algorithm
+        // to check if the current flow is optimal
+        if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
+          typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
+          ShiftCostMap shift_cost(_res_cost, 1);
+          BellmanFord<ResGraph, ShiftCostMap> bf(*_res_graph, shift_cost);
+          bf.init(0);
+          bool done = false;
+          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
+          for (int i = 0; i < K && !done; ++i)
+            done = bf.processNextWeakRound();
+          if (done) break;
+        }
+
+        // Saturate edges not satisfying the optimality condition
+        Capacity delta;
+        for (EdgeIt e(_graph); e != INVALID; ++e) {
+          if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
+            delta = _capacity[e] - (*_flow)[e];
+            _excess[_graph.source(e)] -= delta;
+            _excess[_graph.target(e)] += delta;
+            (*_flow)[e] = _capacity[e];
+          }
+          if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
+            _excess[_graph.target(e)] -= (*_flow)[e];
+            _excess[_graph.source(e)] += (*_flow)[e];
+            (*_flow)[e] = 0;
+          }
+        }
+
+        // Find active nodes (i.e. nodes with positive excess)
+        for (NodeIt n(_graph); n != INVALID; ++n) {
+          if (_excess[n] > 0) active_nodes.push_back(n);
+        }
+
+        // Initialize the next edge maps
+        for (NodeIt n(_graph); n != INVALID; ++n) {
+          next_out[n] = OutEdgeIt(_graph, n);
+          next_in[n] = InEdgeIt(_graph, n);
+          next_dir[n] = true;
+        }
+
+        // Perform push and relabel operations
+        while (active_nodes.size() > 0) {
+          // Select an active node (FIFO selection)
           Node n = active_nodes[0], t;
           bool relabel_enabled = true;
 
-          // Performing push operations if there are admissible edges
-          if (_excess[n] > 0) {
-            for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+          // Perform push operations if there are admissible edges
+          if (_excess[n] > 0 && next_dir[n]) {
+            OutEdgeIt e = next_out[n];
+            for ( ; e != INVALID; ++e) {
               if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
-                delta = _capacity[e] - (*_flow)[e] <= _excess[n] ?
-                        _capacity[e] - (*_flow)[e] : _excess[n];
+                delta = std::min(_capacity[e] - (*_flow)[e], _excess[n]);
                 t = _graph.target(e);
 
@@ -538 +721 @@
                 if (ahead < 0) ahead = 0;
 
-                // Pushing flow along the edge
+                // Push flow along the edge
                 if (ahead < delta) {
                   (*_flow)[e] += ahead;
@@ -558 +741 @@
               }
             }
-          }
-
-          if (_excess[n] > 0) {
-            for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+            if (e != INVALID) {
+              next_out[n] = e;
+            } else {
+              next_dir[n] = false;
+            }
+          }
+
+          if (_excess[n] > 0 && !next_dir[n]) {
+            InEdgeIt e = next_in[n];
+            for ( ; e != INVALID; ++e) {
               if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
-                delta = (*_flow)[e] <= _excess[n] ? (*_flow)[e] : _excess[n];
+                delta = std::min((*_flow)[e], _excess[n]);
                 t = _graph.source(e);
 
@@ -578 +767 @@
                 if (ahead < 0) ahead = 0;
 
-                // Pushing flow along the edge
+                // Push flow along the edge
                 if (ahead < delta) {
                   (*_flow)[e] -= ahead;
@@ -598 +787 @@
               }
             }
-          }
-
+            next_in[n] = e;
+          }
+
+          // Relabel the node if it is still active (or hyper)
           if (relabel_enabled && (_excess[n] > 0 || hyper[n])) {
-            // Performing relabel operation if the node is still active
-            LCost min_red_cost = std::numeric_limits<LCost>::max();
+            LCost min_red_cost = std::numeric_limits<LCost>::max() / 2;
             for (OutEdgeIt oe(_graph, n); oe != INVALID; ++oe) {
               if ( _capacity[oe] - (*_flow)[oe] > 0 &&
@@ -614 +804 @@
             (*_potential)[n] -= min_red_cost + _epsilon;
             hyper[n] = false;
-          }
-
-          // Removing active nodes with non-positive excess
+
+            // Reset the next edge maps
+            next_out[n] = OutEdgeIt(_graph, n);
+            next_in[n] = InEdgeIt(_graph, n);
+            next_dir[n] = true;
+          }
+
+          // Remove nodes that are not active nor hyper
           while ( active_nodes.size() > 0 &&
                   _excess[active_nodes[0]] <= 0 &&
@@ -625 +820 @@
       }
 
-      // Computing node potentials for the original costs
+      // Compute node potentials for the original costs
       ResidualCostMap<CostMap> res_cost(_orig_cost);
       BellmanFord< ResGraph, ResidualCostMap<CostMap> >
@@ -633 +828 @@
         (*_potential)[n] = bf.dist(n);
 
-      // Handling non-zero lower bounds
+      // Handle non-zero lower bounds
       if (_lower) {
         for (EdgeIt e(_graph); e != INVALID; ++e)