# HG changeset patch
# User Peter Kovacs <kpeter@inf.elte.hu>
# Date 2009-11-12 23:45:15
# Node ID fe80a8145653d19e2039c8894311d0f9bd9e704a
# Parent  3b53491bf6436ab241c8ccf8f3c18c529174214a

Small implementation improvements in MCF algorithms (#180)

 - Handle max() as infinite value (not only infinity()).
 - Better GEQ handling in CapacityScaling.
 - Skip the unnecessary saturating operations in the first phase in
   CapacityScaling.
 - Use vector<char> instead of vector<bool> and vector<int> if it is
   possible and it proved to be usually faster.

diff --git a/lemon/capacity_scaling.h b/lemon/capacity_scaling.h
--- a/lemon/capacity_scaling.h
+++ b/lemon/capacity_scaling.h
@@ -133,7 +133,7 @@
     TEMPLATE_DIGRAPH_TYPEDEFS(GR);
 
     typedef std::vector<int> IntVector;
-    typedef std::vector<bool> BoolVector;
+    typedef std::vector<char> BoolVector;
     typedef std::vector<Value> ValueVector;
     typedef std::vector<Cost> CostVector;
 
@@ -196,6 +196,7 @@
     private:
 
       int _node_num;
+      bool _geq;
       const IntVector &_first_out;
       const IntVector &_target;
       const CostVector &_cost;
@@ -210,10 +211,10 @@
     public:
 
       ResidualDijkstra(CapacityScaling& cs) :
-        _node_num(cs._node_num), _first_out(cs._first_out), 
-        _target(cs._target), _cost(cs._cost), _res_cap(cs._res_cap),
-        _excess(cs._excess), _pi(cs._pi), _pred(cs._pred),
-        _dist(cs._node_num)
+        _node_num(cs._node_num), _geq(cs._sum_supply < 0),
+        _first_out(cs._first_out), _target(cs._target), _cost(cs._cost),
+        _res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi),
+        _pred(cs._pred), _dist(cs._node_num)
       {}
 
       int run(int s, Value delta = 1) {
@@ -232,7 +233,8 @@
           heap.pop();
 
           // Traverse outgoing residual arcs
-          for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
+          int last_out = _geq ? _first_out[u+1] : _first_out[u+1] - 1;
+          for (int a = _first_out[u]; a != last_out; ++a) {
             if (_res_cap[a] < delta) continue;
             v = _target[a];
             switch (heap.state(v)) {
@@ -687,22 +689,25 @@
       }
       if (_sum_supply > 0) return INFEASIBLE;
       
-      // Initialize maps
+      // Initialize vectors
       for (int i = 0; i != _root; ++i) {
         _pi[i] = 0;
         _excess[i] = _supply[i];
       }
 
       // Remove non-zero lower bounds
+      const Value MAX = std::numeric_limits<Value>::max();
+      int last_out;
       if (_have_lower) {
         for (int i = 0; i != _root; ++i) {
-          for (int j = _first_out[i]; j != _first_out[i+1]; ++j) {
+          last_out = _first_out[i+1];
+          for (int j = _first_out[i]; j != last_out; ++j) {
             if (_forward[j]) {
               Value c = _lower[j];
               if (c >= 0) {
-                _res_cap[j] = _upper[j] < INF ? _upper[j] - c : INF;
+                _res_cap[j] = _upper[j] < MAX ? _upper[j] - c : INF;
               } else {
-                _res_cap[j] = _upper[j] < INF + c ? _upper[j] - c : INF;
+                _res_cap[j] = _upper[j] < MAX + c ? _upper[j] - c : INF;
               }
               _excess[i] -= c;
               _excess[_target[j]] += c;
@@ -718,15 +723,16 @@
       }
 
       // Handle negative costs
-      for (int u = 0; u != _root; ++u) {
-        for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
-          Value rc = _res_cap[a];
-          if (_cost[a] < 0 && rc > 0) {
-            if (rc == INF) return UNBOUNDED;
-            _excess[u] -= rc;
-            _excess[_target[a]] += rc;
-            _res_cap[a] = 0;
-            _res_cap[_reverse[a]] += rc;
+      for (int i = 0; i != _root; ++i) {
+        last_out = _first_out[i+1] - 1;
+        for (int j = _first_out[i]; j != last_out; ++j) {
+          Value rc = _res_cap[j];
+          if (_cost[j] < 0 && rc > 0) {
+            if (rc >= MAX) return UNBOUNDED;
+            _excess[i] -= rc;
+            _excess[_target[j]] += rc;
+            _res_cap[j] = 0;
+            _res_cap[_reverse[j]] += rc;
           }
         }
       }
@@ -736,24 +742,21 @@
         _pi[_root] = 0;
         _excess[_root] = -_sum_supply;
         for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
-          int u = _target[a];
-          if (_excess[u] < 0) {
-            _res_cap[a] = -_excess[u] + 1;
-          } else {
-            _res_cap[a] = 1;
-          }
-          _res_cap[_reverse[a]] = 0;
+          int ra = _reverse[a];
+          _res_cap[a] = -_sum_supply + 1;
+          _res_cap[ra] = 0;
           _cost[a] = 0;
-          _cost[_reverse[a]] = 0;
+          _cost[ra] = 0;
         }
       } else {
         _pi[_root] = 0;
         _excess[_root] = 0;
         for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
+          int ra = _reverse[a];
           _res_cap[a] = 1;
-          _res_cap[_reverse[a]] = 0;
+          _res_cap[ra] = 0;
           _cost[a] = 0;
-          _cost[_reverse[a]] = 0;
+          _cost[ra] = 0;
         }
       }
 
@@ -762,8 +765,9 @@
         // With scaling
         Value max_sup = 0, max_dem = 0;
         for (int i = 0; i != _node_num; ++i) {
-          if ( _excess[i] > max_sup) max_sup =  _excess[i];
-          if (-_excess[i] > max_dem) max_dem = -_excess[i];
+          Value ex = _excess[i];
+          if ( ex > max_sup) max_sup =  ex;
+          if (-ex > max_dem) max_dem = -ex;
         }
         Value max_cap = 0;
         for (int j = 0; j != _res_arc_num; ++j) {
@@ -789,7 +793,8 @@
 
       // Handle non-zero lower bounds
       if (_have_lower) {
-        for (int j = 0; j != _res_arc_num - _node_num + 1; ++j) {
+        int limit = _first_out[_root];
+        for (int j = 0; j != limit; ++j) {
           if (!_forward[j]) _res_cap[j] += _lower[j];
         }
       }
@@ -812,8 +817,11 @@
       ResidualDijkstra _dijkstra(*this);
       while (true) {
         // Saturate all arcs not satisfying the optimality condition
+        int last_out;
         for (int u = 0; u != _node_num; ++u) {
-          for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
+          last_out = _sum_supply < 0 ?
+            _first_out[u+1] : _first_out[u+1] - 1;
+          for (int a = _first_out[u]; a != last_out; ++a) {
             int v = _target[a];
             Cost c = _cost[a] + _pi[u] - _pi[v];
             Value rc = _res_cap[a];
@@ -830,8 +838,9 @@
         _excess_nodes.clear();
         _deficit_nodes.clear();
         for (int u = 0; u != _node_num; ++u) {
-          if (_excess[u] >=  _delta) _excess_nodes.push_back(u);
-          if (_excess[u] <= -_delta) _deficit_nodes.push_back(u);
+          Value ex = _excess[u];
+          if (ex >=  _delta) _excess_nodes.push_back(u);
+          if (ex <= -_delta) _deficit_nodes.push_back(u);
         }
         int next_node = 0, next_def_node = 0;
 
diff --git a/lemon/network_simplex.h b/lemon/network_simplex.h
--- a/lemon/network_simplex.h
+++ b/lemon/network_simplex.h
@@ -164,7 +164,7 @@
     TEMPLATE_DIGRAPH_TYPEDEFS(GR);
 
     typedef std::vector<int> IntVector;
-    typedef std::vector<bool> BoolVector;
+    typedef std::vector<char> CharVector;
     typedef std::vector<Value> ValueVector;
     typedef std::vector<Cost> CostVector;
 
@@ -212,8 +212,8 @@
     IntVector _succ_num;
     IntVector _last_succ;
     IntVector _dirty_revs;
-    BoolVector _forward;
-    IntVector _state;
+    CharVector _forward;
+    CharVector _state;
     int _root;
 
     // Temporary data used in the current pivot iteration
@@ -221,6 +221,8 @@
     int first, second, right, last;
     int stem, par_stem, new_stem;
     Value delta;
+    
+    const Value MAX;
 
   public:
   
@@ -242,7 +244,7 @@
       const IntVector  &_source;
       const IntVector  &_target;
       const CostVector &_cost;
-      const IntVector  &_state;
+      const CharVector &_state;
       const CostVector &_pi;
       int &_in_arc;
       int _search_arc_num;
@@ -294,7 +296,7 @@
       const IntVector  &_source;
       const IntVector  &_target;
       const CostVector &_cost;
-      const IntVector  &_state;
+      const CharVector &_state;
       const CostVector &_pi;
       int &_in_arc;
       int _search_arc_num;
@@ -333,7 +335,7 @@
       const IntVector  &_source;
       const IntVector  &_target;
       const CostVector &_cost;
-      const IntVector  &_state;
+      const CharVector &_state;
       const CostVector &_pi;
       int &_in_arc;
       int _search_arc_num;
@@ -406,7 +408,7 @@
       const IntVector  &_source;
       const IntVector  &_target;
       const CostVector &_cost;
-      const IntVector  &_state;
+      const CharVector &_state;
       const CostVector &_pi;
       int &_in_arc;
       int _search_arc_num;
@@ -509,7 +511,7 @@
       const IntVector  &_source;
       const IntVector  &_target;
       const CostVector &_cost;
-      const IntVector  &_state;
+      const CharVector &_state;
       const CostVector &_pi;
       int &_in_arc;
       int _search_arc_num;
@@ -631,9 +633,9 @@
     /// but it is usually slower. Therefore it is disabled by default.
     NetworkSimplex(const GR& graph, bool arc_mixing = false) :
       _graph(graph), _node_id(graph), _arc_id(graph),
+      MAX(std::numeric_limits<Value>::max()),
       INF(std::numeric_limits<Value>::has_infinity ?
-          std::numeric_limits<Value>::infinity() :
-          std::numeric_limits<Value>::max())
+          std::numeric_limits<Value>::infinity() : MAX)
     {
       // Check the value types
       LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
@@ -1020,9 +1022,9 @@
         for (int i = 0; i != _arc_num; ++i) {
           Value c = _lower[i];
           if (c >= 0) {
-            _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
+            _cap[i] = _upper[i] < MAX ? _upper[i] - c : INF;
           } else {
-            _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
+            _cap[i] = _upper[i] < MAX + c ? _upper[i] - c : INF;
           }
           _supply[_source[i]] -= c;
           _supply[_target[i]] += c;
@@ -1214,7 +1216,7 @@
       for (int u = first; u != join; u = _parent[u]) {
         e = _pred[u];
         d = _forward[u] ?
-          _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
+          _flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]);
         if (d < delta) {
           delta = d;
           u_out = u;
@@ -1225,7 +1227,7 @@
       for (int u = second; u != join; u = _parent[u]) {
         e = _pred[u];
         d = _forward[u] ? 
-          (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
+          (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
         if (d <= delta) {
           delta = d;
           u_out = u;
@@ -1424,7 +1426,7 @@
       while (pivot.findEnteringArc()) {
         findJoinNode();
         bool change = findLeavingArc();
-        if (delta >= INF) return UNBOUNDED;
+        if (delta >= MAX) return UNBOUNDED;
         changeFlow(change);
         if (change) {
           updateTreeStructure();