Small implementation improvements in MCF algorithms (#180)
authorPeter Kovacs <kpeter@inf.elte.hu>
Thu, 12 Nov 2009 23:45:15 +0100
changeset 877fe80a8145653
parent 876 3b53491bf643
child 878 4b1b378823dc
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.
lemon/capacity_scaling.h
lemon/network_simplex.h
     1.1 --- a/lemon/capacity_scaling.h	Thu Nov 12 23:34:35 2009 +0100
     1.2 +++ b/lemon/capacity_scaling.h	Thu Nov 12 23:45:15 2009 +0100
     1.3 @@ -133,7 +133,7 @@
     1.4      TEMPLATE_DIGRAPH_TYPEDEFS(GR);
     1.5  
     1.6      typedef std::vector<int> IntVector;
     1.7 -    typedef std::vector<bool> BoolVector;
     1.8 +    typedef std::vector<char> BoolVector;
     1.9      typedef std::vector<Value> ValueVector;
    1.10      typedef std::vector<Cost> CostVector;
    1.11  
    1.12 @@ -196,6 +196,7 @@
    1.13      private:
    1.14  
    1.15        int _node_num;
    1.16 +      bool _geq;
    1.17        const IntVector &_first_out;
    1.18        const IntVector &_target;
    1.19        const CostVector &_cost;
    1.20 @@ -210,10 +211,10 @@
    1.21      public:
    1.22  
    1.23        ResidualDijkstra(CapacityScaling& cs) :
    1.24 -        _node_num(cs._node_num), _first_out(cs._first_out), 
    1.25 -        _target(cs._target), _cost(cs._cost), _res_cap(cs._res_cap),
    1.26 -        _excess(cs._excess), _pi(cs._pi), _pred(cs._pred),
    1.27 -        _dist(cs._node_num)
    1.28 +        _node_num(cs._node_num), _geq(cs._sum_supply < 0),
    1.29 +        _first_out(cs._first_out), _target(cs._target), _cost(cs._cost),
    1.30 +        _res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi),
    1.31 +        _pred(cs._pred), _dist(cs._node_num)
    1.32        {}
    1.33  
    1.34        int run(int s, Value delta = 1) {
    1.35 @@ -232,7 +233,8 @@
    1.36            heap.pop();
    1.37  
    1.38            // Traverse outgoing residual arcs
    1.39 -          for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
    1.40 +          int last_out = _geq ? _first_out[u+1] : _first_out[u+1] - 1;
    1.41 +          for (int a = _first_out[u]; a != last_out; ++a) {
    1.42              if (_res_cap[a] < delta) continue;
    1.43              v = _target[a];
    1.44              switch (heap.state(v)) {
    1.45 @@ -687,22 +689,25 @@
    1.46        }
    1.47        if (_sum_supply > 0) return INFEASIBLE;
    1.48        
    1.49 -      // Initialize maps
    1.50 +      // Initialize vectors
    1.51        for (int i = 0; i != _root; ++i) {
    1.52          _pi[i] = 0;
    1.53          _excess[i] = _supply[i];
    1.54        }
    1.55  
    1.56        // Remove non-zero lower bounds
    1.57 +      const Value MAX = std::numeric_limits<Value>::max();
    1.58 +      int last_out;
    1.59        if (_have_lower) {
    1.60          for (int i = 0; i != _root; ++i) {
    1.61 -          for (int j = _first_out[i]; j != _first_out[i+1]; ++j) {
    1.62 +          last_out = _first_out[i+1];
    1.63 +          for (int j = _first_out[i]; j != last_out; ++j) {
    1.64              if (_forward[j]) {
    1.65                Value c = _lower[j];
    1.66                if (c >= 0) {
    1.67 -                _res_cap[j] = _upper[j] < INF ? _upper[j] - c : INF;
    1.68 +                _res_cap[j] = _upper[j] < MAX ? _upper[j] - c : INF;
    1.69                } else {
    1.70 -                _res_cap[j] = _upper[j] < INF + c ? _upper[j] - c : INF;
    1.71 +                _res_cap[j] = _upper[j] < MAX + c ? _upper[j] - c : INF;
    1.72                }
    1.73                _excess[i] -= c;
    1.74                _excess[_target[j]] += c;
    1.75 @@ -718,15 +723,16 @@
    1.76        }
    1.77  
    1.78        // Handle negative costs
    1.79 -      for (int u = 0; u != _root; ++u) {
    1.80 -        for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
    1.81 -          Value rc = _res_cap[a];
    1.82 -          if (_cost[a] < 0 && rc > 0) {
    1.83 -            if (rc == INF) return UNBOUNDED;
    1.84 -            _excess[u] -= rc;
    1.85 -            _excess[_target[a]] += rc;
    1.86 -            _res_cap[a] = 0;
    1.87 -            _res_cap[_reverse[a]] += rc;
    1.88 +      for (int i = 0; i != _root; ++i) {
    1.89 +        last_out = _first_out[i+1] - 1;
    1.90 +        for (int j = _first_out[i]; j != last_out; ++j) {
    1.91 +          Value rc = _res_cap[j];
    1.92 +          if (_cost[j] < 0 && rc > 0) {
    1.93 +            if (rc >= MAX) return UNBOUNDED;
    1.94 +            _excess[i] -= rc;
    1.95 +            _excess[_target[j]] += rc;
    1.96 +            _res_cap[j] = 0;
    1.97 +            _res_cap[_reverse[j]] += rc;
    1.98            }
    1.99          }
   1.100        }
   1.101 @@ -736,24 +742,21 @@
   1.102          _pi[_root] = 0;
   1.103          _excess[_root] = -_sum_supply;
   1.104          for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
   1.105 -          int u = _target[a];
   1.106 -          if (_excess[u] < 0) {
   1.107 -            _res_cap[a] = -_excess[u] + 1;
   1.108 -          } else {
   1.109 -            _res_cap[a] = 1;
   1.110 -          }
   1.111 -          _res_cap[_reverse[a]] = 0;
   1.112 +          int ra = _reverse[a];
   1.113 +          _res_cap[a] = -_sum_supply + 1;
   1.114 +          _res_cap[ra] = 0;
   1.115            _cost[a] = 0;
   1.116 -          _cost[_reverse[a]] = 0;
   1.117 +          _cost[ra] = 0;
   1.118          }
   1.119        } else {
   1.120          _pi[_root] = 0;
   1.121          _excess[_root] = 0;
   1.122          for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
   1.123 +          int ra = _reverse[a];
   1.124            _res_cap[a] = 1;
   1.125 -          _res_cap[_reverse[a]] = 0;
   1.126 +          _res_cap[ra] = 0;
   1.127            _cost[a] = 0;
   1.128 -          _cost[_reverse[a]] = 0;
   1.129 +          _cost[ra] = 0;
   1.130          }
   1.131        }
   1.132  
   1.133 @@ -762,8 +765,9 @@
   1.134          // With scaling
   1.135          Value max_sup = 0, max_dem = 0;
   1.136          for (int i = 0; i != _node_num; ++i) {
   1.137 -          if ( _excess[i] > max_sup) max_sup =  _excess[i];
   1.138 -          if (-_excess[i] > max_dem) max_dem = -_excess[i];
   1.139 +          Value ex = _excess[i];
   1.140 +          if ( ex > max_sup) max_sup =  ex;
   1.141 +          if (-ex > max_dem) max_dem = -ex;
   1.142          }
   1.143          Value max_cap = 0;
   1.144          for (int j = 0; j != _res_arc_num; ++j) {
   1.145 @@ -789,7 +793,8 @@
   1.146  
   1.147        // Handle non-zero lower bounds
   1.148        if (_have_lower) {
   1.149 -        for (int j = 0; j != _res_arc_num - _node_num + 1; ++j) {
   1.150 +        int limit = _first_out[_root];
   1.151 +        for (int j = 0; j != limit; ++j) {
   1.152            if (!_forward[j]) _res_cap[j] += _lower[j];
   1.153          }
   1.154        }
   1.155 @@ -812,8 +817,11 @@
   1.156        ResidualDijkstra _dijkstra(*this);
   1.157        while (true) {
   1.158          // Saturate all arcs not satisfying the optimality condition
   1.159 +        int last_out;
   1.160          for (int u = 0; u != _node_num; ++u) {
   1.161 -          for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
   1.162 +          last_out = _sum_supply < 0 ?
   1.163 +            _first_out[u+1] : _first_out[u+1] - 1;
   1.164 +          for (int a = _first_out[u]; a != last_out; ++a) {
   1.165              int v = _target[a];
   1.166              Cost c = _cost[a] + _pi[u] - _pi[v];
   1.167              Value rc = _res_cap[a];
   1.168 @@ -830,8 +838,9 @@
   1.169          _excess_nodes.clear();
   1.170          _deficit_nodes.clear();
   1.171          for (int u = 0; u != _node_num; ++u) {
   1.172 -          if (_excess[u] >=  _delta) _excess_nodes.push_back(u);
   1.173 -          if (_excess[u] <= -_delta) _deficit_nodes.push_back(u);
   1.174 +          Value ex = _excess[u];
   1.175 +          if (ex >=  _delta) _excess_nodes.push_back(u);
   1.176 +          if (ex <= -_delta) _deficit_nodes.push_back(u);
   1.177          }
   1.178          int next_node = 0, next_def_node = 0;
   1.179  
     2.1 --- a/lemon/network_simplex.h	Thu Nov 12 23:34:35 2009 +0100
     2.2 +++ b/lemon/network_simplex.h	Thu Nov 12 23:45:15 2009 +0100
     2.3 @@ -164,7 +164,7 @@
     2.4      TEMPLATE_DIGRAPH_TYPEDEFS(GR);
     2.5  
     2.6      typedef std::vector<int> IntVector;
     2.7 -    typedef std::vector<bool> BoolVector;
     2.8 +    typedef std::vector<char> CharVector;
     2.9      typedef std::vector<Value> ValueVector;
    2.10      typedef std::vector<Cost> CostVector;
    2.11  
    2.12 @@ -212,8 +212,8 @@
    2.13      IntVector _succ_num;
    2.14      IntVector _last_succ;
    2.15      IntVector _dirty_revs;
    2.16 -    BoolVector _forward;
    2.17 -    IntVector _state;
    2.18 +    CharVector _forward;
    2.19 +    CharVector _state;
    2.20      int _root;
    2.21  
    2.22      // Temporary data used in the current pivot iteration
    2.23 @@ -221,6 +221,8 @@
    2.24      int first, second, right, last;
    2.25      int stem, par_stem, new_stem;
    2.26      Value delta;
    2.27 +    
    2.28 +    const Value MAX;
    2.29  
    2.30    public:
    2.31    
    2.32 @@ -242,7 +244,7 @@
    2.33        const IntVector  &_source;
    2.34        const IntVector  &_target;
    2.35        const CostVector &_cost;
    2.36 -      const IntVector  &_state;
    2.37 +      const CharVector &_state;
    2.38        const CostVector &_pi;
    2.39        int &_in_arc;
    2.40        int _search_arc_num;
    2.41 @@ -294,7 +296,7 @@
    2.42        const IntVector  &_source;
    2.43        const IntVector  &_target;
    2.44        const CostVector &_cost;
    2.45 -      const IntVector  &_state;
    2.46 +      const CharVector &_state;
    2.47        const CostVector &_pi;
    2.48        int &_in_arc;
    2.49        int _search_arc_num;
    2.50 @@ -333,7 +335,7 @@
    2.51        const IntVector  &_source;
    2.52        const IntVector  &_target;
    2.53        const CostVector &_cost;
    2.54 -      const IntVector  &_state;
    2.55 +      const CharVector &_state;
    2.56        const CostVector &_pi;
    2.57        int &_in_arc;
    2.58        int _search_arc_num;
    2.59 @@ -406,7 +408,7 @@
    2.60        const IntVector  &_source;
    2.61        const IntVector  &_target;
    2.62        const CostVector &_cost;
    2.63 -      const IntVector  &_state;
    2.64 +      const CharVector &_state;
    2.65        const CostVector &_pi;
    2.66        int &_in_arc;
    2.67        int _search_arc_num;
    2.68 @@ -509,7 +511,7 @@
    2.69        const IntVector  &_source;
    2.70        const IntVector  &_target;
    2.71        const CostVector &_cost;
    2.72 -      const IntVector  &_state;
    2.73 +      const CharVector &_state;
    2.74        const CostVector &_pi;
    2.75        int &_in_arc;
    2.76        int _search_arc_num;
    2.77 @@ -631,9 +633,9 @@
    2.78      /// but it is usually slower. Therefore it is disabled by default.
    2.79      NetworkSimplex(const GR& graph, bool arc_mixing = false) :
    2.80        _graph(graph), _node_id(graph), _arc_id(graph),
    2.81 +      MAX(std::numeric_limits<Value>::max()),
    2.82        INF(std::numeric_limits<Value>::has_infinity ?
    2.83 -          std::numeric_limits<Value>::infinity() :
    2.84 -          std::numeric_limits<Value>::max())
    2.85 +          std::numeric_limits<Value>::infinity() : MAX)
    2.86      {
    2.87        // Check the value types
    2.88        LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
    2.89 @@ -1020,9 +1022,9 @@
    2.90          for (int i = 0; i != _arc_num; ++i) {
    2.91            Value c = _lower[i];
    2.92            if (c >= 0) {
    2.93 -            _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
    2.94 +            _cap[i] = _upper[i] < MAX ? _upper[i] - c : INF;
    2.95            } else {
    2.96 -            _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
    2.97 +            _cap[i] = _upper[i] < MAX + c ? _upper[i] - c : INF;
    2.98            }
    2.99            _supply[_source[i]] -= c;
   2.100            _supply[_target[i]] += c;
   2.101 @@ -1214,7 +1216,7 @@
   2.102        for (int u = first; u != join; u = _parent[u]) {
   2.103          e = _pred[u];
   2.104          d = _forward[u] ?
   2.105 -          _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
   2.106 +          _flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]);
   2.107          if (d < delta) {
   2.108            delta = d;
   2.109            u_out = u;
   2.110 @@ -1225,7 +1227,7 @@
   2.111        for (int u = second; u != join; u = _parent[u]) {
   2.112          e = _pred[u];
   2.113          d = _forward[u] ? 
   2.114 -          (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
   2.115 +          (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
   2.116          if (d <= delta) {
   2.117            delta = d;
   2.118            u_out = u;
   2.119 @@ -1424,7 +1426,7 @@
   2.120        while (pivot.findEnteringArc()) {
   2.121          findJoinNode();
   2.122          bool change = findLeavingArc();
   2.123 -        if (delta >= INF) return UNBOUNDED;
   2.124 +        if (delta >= MAX) return UNBOUNDED;
   2.125          changeFlow(change);
   2.126          if (change) {
   2.127            updateTreeStructure();