Merge #340
authorAlpar Juttner <alpar@cs.elte.hu>
Fri, 26 Feb 2010 14:00:20 +0100
changeset 8402914b6f0fde0
parent 838 2c35bef44dd1
parent 839 f3bc4e9b5f3a
child 841 aa8c9008b3de
Merge #340
lemon/capacity_scaling.h
lemon/cost_scaling.h
lemon/cycle_canceling.h
lemon/network_simplex.h
     1.1 --- a/lemon/capacity_scaling.h	Sun Feb 21 18:55:30 2010 +0100
     1.2 +++ b/lemon/capacity_scaling.h	Fri Feb 26 14:00:20 2010 +0100
     1.3 @@ -139,9 +139,10 @@
     1.4      TEMPLATE_DIGRAPH_TYPEDEFS(GR);
     1.5  
     1.6      typedef std::vector<int> IntVector;
     1.7 -    typedef std::vector<char> BoolVector;
     1.8      typedef std::vector<Value> ValueVector;
     1.9      typedef std::vector<Cost> CostVector;
    1.10 +    typedef std::vector<char> BoolVector;
    1.11 +    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
    1.12  
    1.13    private:
    1.14  
    1.15 @@ -798,15 +799,15 @@
    1.16        // Initialize delta value
    1.17        if (_factor > 1) {
    1.18          // With scaling
    1.19 -        Value max_sup = 0, max_dem = 0;
    1.20 -        for (int i = 0; i != _node_num; ++i) {
    1.21 +        Value max_sup = 0, max_dem = 0, max_cap = 0;
    1.22 +        for (int i = 0; i != _root; ++i) {
    1.23            Value ex = _excess[i];
    1.24            if ( ex > max_sup) max_sup =  ex;
    1.25            if (-ex > max_dem) max_dem = -ex;
    1.26 -        }
    1.27 -        Value max_cap = 0;
    1.28 -        for (int j = 0; j != _res_arc_num; ++j) {
    1.29 -          if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
    1.30 +          int last_out = _first_out[i+1] - 1;
    1.31 +          for (int j = _first_out[i]; j != last_out; ++j) {
    1.32 +            if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
    1.33 +          }
    1.34          }
    1.35          max_sup = std::min(std::min(max_sup, max_dem), max_cap);
    1.36          for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) ;
     2.1 --- a/lemon/cost_scaling.h	Sun Feb 21 18:55:30 2010 +0100
     2.2 +++ b/lemon/cost_scaling.h	Fri Feb 26 14:00:20 2010 +0100
     2.3 @@ -201,10 +201,11 @@
     2.4      TEMPLATE_DIGRAPH_TYPEDEFS(GR);
     2.5  
     2.6      typedef std::vector<int> IntVector;
     2.7 -    typedef std::vector<char> BoolVector;
     2.8      typedef std::vector<Value> ValueVector;
     2.9      typedef std::vector<Cost> CostVector;
    2.10      typedef std::vector<LargeCost> LargeCostVector;
    2.11 +    typedef std::vector<char> BoolVector;
    2.12 +    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
    2.13  
    2.14    private:
    2.15    
    2.16 @@ -248,6 +249,7 @@
    2.17      // Parameters of the problem
    2.18      bool _have_lower;
    2.19      Value _sum_supply;
    2.20 +    int _sup_node_num;
    2.21  
    2.22      // Data structures for storing the digraph
    2.23      IntNodeMap _node_id;
    2.24 @@ -276,6 +278,12 @@
    2.25      LargeCost _epsilon;
    2.26      int _alpha;
    2.27  
    2.28 +    IntVector _buckets;
    2.29 +    IntVector _bucket_next;
    2.30 +    IntVector _bucket_prev;
    2.31 +    IntVector _rank;
    2.32 +    int _max_rank;
    2.33 +  
    2.34      // Data for a StaticDigraph structure
    2.35      typedef std::pair<int, int> IntPair;
    2.36      StaticDigraph _sgr;
    2.37 @@ -828,6 +836,11 @@
    2.38          }
    2.39        }
    2.40  
    2.41 +      _sup_node_num = 0;
    2.42 +      for (NodeIt n(_graph); n != INVALID; ++n) {
    2.43 +        if (sup[n] > 0) ++_sup_node_num;
    2.44 +      }
    2.45 +
    2.46        // Find a feasible flow using Circulation
    2.47        Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
    2.48          circ(_graph, low, cap, sup);
    2.49 @@ -862,7 +875,7 @@
    2.50          }
    2.51          for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
    2.52            int ra = _reverse[a];
    2.53 -          _res_cap[a] = 1;
    2.54 +          _res_cap[a] = 0;
    2.55            _res_cap[ra] = 0;
    2.56            _cost[a] = 0;
    2.57            _cost[ra] = 0;
    2.58 @@ -876,7 +889,14 @@
    2.59      void start(Method method) {
    2.60        // Maximum path length for partial augment
    2.61        const int MAX_PATH_LENGTH = 4;
    2.62 -      
    2.63 +
    2.64 +      // Initialize data structures for buckets      
    2.65 +      _max_rank = _alpha * _res_node_num;
    2.66 +      _buckets.resize(_max_rank);
    2.67 +      _bucket_next.resize(_res_node_num + 1);
    2.68 +      _bucket_prev.resize(_res_node_num + 1);
    2.69 +      _rank.resize(_res_node_num + 1);
    2.70 +  
    2.71        // Execute the algorithm
    2.72        switch (method) {
    2.73          case PUSH:
    2.74 @@ -915,63 +935,175 @@
    2.75          }
    2.76        }
    2.77      }
    2.78 +    
    2.79 +    // Initialize a cost scaling phase
    2.80 +    void initPhase() {
    2.81 +      // Saturate arcs not satisfying the optimality condition
    2.82 +      for (int u = 0; u != _res_node_num; ++u) {
    2.83 +        int last_out = _first_out[u+1];
    2.84 +        LargeCost pi_u = _pi[u];
    2.85 +        for (int a = _first_out[u]; a != last_out; ++a) {
    2.86 +          int v = _target[a];
    2.87 +          if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) {
    2.88 +            Value delta = _res_cap[a];
    2.89 +            _excess[u] -= delta;
    2.90 +            _excess[v] += delta;
    2.91 +            _res_cap[a] = 0;
    2.92 +            _res_cap[_reverse[a]] += delta;
    2.93 +          }
    2.94 +        }
    2.95 +      }
    2.96 +      
    2.97 +      // Find active nodes (i.e. nodes with positive excess)
    2.98 +      for (int u = 0; u != _res_node_num; ++u) {
    2.99 +        if (_excess[u] > 0) _active_nodes.push_back(u);
   2.100 +      }
   2.101 +
   2.102 +      // Initialize the next arcs
   2.103 +      for (int u = 0; u != _res_node_num; ++u) {
   2.104 +        _next_out[u] = _first_out[u];
   2.105 +      }
   2.106 +    }
   2.107 +    
   2.108 +    // Early termination heuristic
   2.109 +    bool earlyTermination() {
   2.110 +      const double EARLY_TERM_FACTOR = 3.0;
   2.111 +
   2.112 +      // Build a static residual graph
   2.113 +      _arc_vec.clear();
   2.114 +      _cost_vec.clear();
   2.115 +      for (int j = 0; j != _res_arc_num; ++j) {
   2.116 +        if (_res_cap[j] > 0) {
   2.117 +          _arc_vec.push_back(IntPair(_source[j], _target[j]));
   2.118 +          _cost_vec.push_back(_cost[j] + 1);
   2.119 +        }
   2.120 +      }
   2.121 +      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
   2.122 +
   2.123 +      // Run Bellman-Ford algorithm to check if the current flow is optimal
   2.124 +      BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
   2.125 +      bf.init(0);
   2.126 +      bool done = false;
   2.127 +      int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num)));
   2.128 +      for (int i = 0; i < K && !done; ++i) {
   2.129 +        done = bf.processNextWeakRound();
   2.130 +      }
   2.131 +      return done;
   2.132 +    }
   2.133 +
   2.134 +    // Global potential update heuristic
   2.135 +    void globalUpdate() {
   2.136 +      int bucket_end = _root + 1;
   2.137 +    
   2.138 +      // Initialize buckets
   2.139 +      for (int r = 0; r != _max_rank; ++r) {
   2.140 +        _buckets[r] = bucket_end;
   2.141 +      }
   2.142 +      Value total_excess = 0;
   2.143 +      for (int i = 0; i != _res_node_num; ++i) {
   2.144 +        if (_excess[i] < 0) {
   2.145 +          _rank[i] = 0;
   2.146 +          _bucket_next[i] = _buckets[0];
   2.147 +          _bucket_prev[_buckets[0]] = i;
   2.148 +          _buckets[0] = i;
   2.149 +        } else {
   2.150 +          total_excess += _excess[i];
   2.151 +          _rank[i] = _max_rank;
   2.152 +        }
   2.153 +      }
   2.154 +      if (total_excess == 0) return;
   2.155 +
   2.156 +      // Search the buckets
   2.157 +      int r = 0;
   2.158 +      for ( ; r != _max_rank; ++r) {
   2.159 +        while (_buckets[r] != bucket_end) {
   2.160 +          // Remove the first node from the current bucket
   2.161 +          int u = _buckets[r];
   2.162 +          _buckets[r] = _bucket_next[u];
   2.163 +          
   2.164 +          // Search the incomming arcs of u
   2.165 +          LargeCost pi_u = _pi[u];
   2.166 +          int last_out = _first_out[u+1];
   2.167 +          for (int a = _first_out[u]; a != last_out; ++a) {
   2.168 +            int ra = _reverse[a];
   2.169 +            if (_res_cap[ra] > 0) {
   2.170 +              int v = _source[ra];
   2.171 +              int old_rank_v = _rank[v];
   2.172 +              if (r < old_rank_v) {
   2.173 +                // Compute the new rank of v
   2.174 +                LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
   2.175 +                int new_rank_v = old_rank_v;
   2.176 +                if (nrc < LargeCost(_max_rank))
   2.177 +                  new_rank_v = r + 1 + int(nrc);
   2.178 +                  
   2.179 +                // Change the rank of v
   2.180 +                if (new_rank_v < old_rank_v) {
   2.181 +                  _rank[v] = new_rank_v;
   2.182 +                  _next_out[v] = _first_out[v];
   2.183 +                  
   2.184 +                  // Remove v from its old bucket
   2.185 +                  if (old_rank_v < _max_rank) {
   2.186 +                    if (_buckets[old_rank_v] == v) {
   2.187 +                      _buckets[old_rank_v] = _bucket_next[v];
   2.188 +                    } else {
   2.189 +                      _bucket_next[_bucket_prev[v]] = _bucket_next[v];
   2.190 +                      _bucket_prev[_bucket_next[v]] = _bucket_prev[v];
   2.191 +                    }
   2.192 +                  }
   2.193 +                  
   2.194 +                  // Insert v to its new bucket
   2.195 +                  _bucket_next[v] = _buckets[new_rank_v];
   2.196 +                  _bucket_prev[_buckets[new_rank_v]] = v;
   2.197 +                  _buckets[new_rank_v] = v;
   2.198 +                }
   2.199 +              }
   2.200 +            }
   2.201 +          }
   2.202 +
   2.203 +          // Finish search if there are no more active nodes
   2.204 +          if (_excess[u] > 0) {
   2.205 +            total_excess -= _excess[u];
   2.206 +            if (total_excess <= 0) break;
   2.207 +          }
   2.208 +        }
   2.209 +        if (total_excess <= 0) break;
   2.210 +      }
   2.211 +      
   2.212 +      // Relabel nodes
   2.213 +      for (int u = 0; u != _res_node_num; ++u) {
   2.214 +        int k = std::min(_rank[u], r);
   2.215 +        if (k > 0) {
   2.216 +          _pi[u] -= _epsilon * k;
   2.217 +          _next_out[u] = _first_out[u];
   2.218 +        }
   2.219 +      }
   2.220 +    }
   2.221  
   2.222      /// Execute the algorithm performing augment and relabel operations
   2.223      void startAugment(int max_length = std::numeric_limits<int>::max()) {
   2.224        // Paramters for heuristics
   2.225 -      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
   2.226 -      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
   2.227 +      const int EARLY_TERM_EPSILON_LIMIT = 1000;
   2.228 +      const double GLOBAL_UPDATE_FACTOR = 3.0;
   2.229  
   2.230 +      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
   2.231 +        (_res_node_num + _sup_node_num * _sup_node_num));
   2.232 +      int next_update_limit = global_update_freq;
   2.233 +      
   2.234 +      int relabel_cnt = 0;
   2.235 +      
   2.236        // Perform cost scaling phases
   2.237 -      IntVector pred_arc(_res_node_num);
   2.238 -      std::vector<int> path_nodes;
   2.239 +      std::vector<int> path;
   2.240        for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
   2.241                                          1 : _epsilon / _alpha )
   2.242        {
   2.243 -        // "Early Termination" heuristic: use Bellman-Ford algorithm
   2.244 -        // to check if the current flow is optimal
   2.245 -        if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
   2.246 -          _arc_vec.clear();
   2.247 -          _cost_vec.clear();
   2.248 -          for (int j = 0; j != _res_arc_num; ++j) {
   2.249 -            if (_res_cap[j] > 0) {
   2.250 -              _arc_vec.push_back(IntPair(_source[j], _target[j]));
   2.251 -              _cost_vec.push_back(_cost[j] + 1);
   2.252 -            }
   2.253 -          }
   2.254 -          _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
   2.255 -
   2.256 -          BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
   2.257 -          bf.init(0);
   2.258 -          bool done = false;
   2.259 -          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
   2.260 -          for (int i = 0; i < K && !done; ++i)
   2.261 -            done = bf.processNextWeakRound();
   2.262 -          if (done) break;
   2.263 -        }
   2.264 -
   2.265 -        // Saturate arcs not satisfying the optimality condition
   2.266 -        for (int a = 0; a != _res_arc_num; ++a) {
   2.267 -          if (_res_cap[a] > 0 &&
   2.268 -              _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
   2.269 -            Value delta = _res_cap[a];
   2.270 -            _excess[_source[a]] -= delta;
   2.271 -            _excess[_target[a]] += delta;
   2.272 -            _res_cap[a] = 0;
   2.273 -            _res_cap[_reverse[a]] += delta;
   2.274 -          }
   2.275 +        // Early termination heuristic
   2.276 +        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
   2.277 +          if (earlyTermination()) break;
   2.278          }
   2.279          
   2.280 -        // Find active nodes (i.e. nodes with positive excess)
   2.281 -        for (int u = 0; u != _res_node_num; ++u) {
   2.282 -          if (_excess[u] > 0) _active_nodes.push_back(u);
   2.283 -        }
   2.284 -
   2.285 -        // Initialize the next arcs
   2.286 -        for (int u = 0; u != _res_node_num; ++u) {
   2.287 -          _next_out[u] = _first_out[u];
   2.288 -        }
   2.289 -
   2.290 +        // Initialize current phase
   2.291 +        initPhase();
   2.292 +        
   2.293          // Perform partial augment and relabel operations
   2.294          while (true) {
   2.295            // Select an active node (FIFO selection)
   2.296 @@ -981,46 +1113,44 @@
   2.297            }
   2.298            if (_active_nodes.size() == 0) break;
   2.299            int start = _active_nodes.front();
   2.300 -          path_nodes.clear();
   2.301 -          path_nodes.push_back(start);
   2.302  
   2.303            // Find an augmenting path from the start node
   2.304 +          path.clear();
   2.305            int tip = start;
   2.306 -          while (_excess[tip] >= 0 &&
   2.307 -                 int(path_nodes.size()) <= max_length) {
   2.308 +          while (_excess[tip] >= 0 && int(path.size()) < max_length) {
   2.309              int u;
   2.310 -            LargeCost min_red_cost, rc;
   2.311 -            int last_out = _sum_supply < 0 ?
   2.312 -              _first_out[tip+1] : _first_out[tip+1] - 1;
   2.313 +            LargeCost min_red_cost, rc, pi_tip = _pi[tip];
   2.314 +            int last_out = _first_out[tip+1];
   2.315              for (int a = _next_out[tip]; a != last_out; ++a) {
   2.316 -              if (_res_cap[a] > 0 &&
   2.317 -                  _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
   2.318 -                u = _target[a];
   2.319 -                pred_arc[u] = a;
   2.320 +              u = _target[a];
   2.321 +              if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) {
   2.322 +                path.push_back(a);
   2.323                  _next_out[tip] = a;
   2.324                  tip = u;
   2.325 -                path_nodes.push_back(tip);
   2.326                  goto next_step;
   2.327                }
   2.328              }
   2.329  
   2.330              // Relabel tip node
   2.331 -            min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
   2.332 +            min_red_cost = std::numeric_limits<LargeCost>::max();
   2.333 +            if (tip != start) {
   2.334 +              int ra = _reverse[path.back()];
   2.335 +              min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]];
   2.336 +            }
   2.337              for (int a = _first_out[tip]; a != last_out; ++a) {
   2.338 -              rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
   2.339 +              rc = _cost[a] + pi_tip - _pi[_target[a]];
   2.340                if (_res_cap[a] > 0 && rc < min_red_cost) {
   2.341                  min_red_cost = rc;
   2.342                }
   2.343              }
   2.344              _pi[tip] -= min_red_cost + _epsilon;
   2.345 -
   2.346 -            // Reset the next arc of tip
   2.347              _next_out[tip] = _first_out[tip];
   2.348 +            ++relabel_cnt;
   2.349  
   2.350              // Step back
   2.351              if (tip != start) {
   2.352 -              path_nodes.pop_back();
   2.353 -              tip = path_nodes.back();
   2.354 +              tip = _source[path.back()];
   2.355 +              path.pop_back();
   2.356              }
   2.357  
   2.358            next_step: ;
   2.359 @@ -1028,11 +1158,11 @@
   2.360  
   2.361            // Augment along the found path (as much flow as possible)
   2.362            Value delta;
   2.363 -          int u, v = path_nodes.front(), pa;
   2.364 -          for (int i = 1; i < int(path_nodes.size()); ++i) {
   2.365 +          int pa, u, v = start;
   2.366 +          for (int i = 0; i != int(path.size()); ++i) {
   2.367 +            pa = path[i];
   2.368              u = v;
   2.369 -            v = path_nodes[i];
   2.370 -            pa = pred_arc[v];
   2.371 +            v = _target[pa];
   2.372              delta = std::min(_res_cap[pa], _excess[u]);
   2.373              _res_cap[pa] -= delta;
   2.374              _res_cap[_reverse[pa]] += delta;
   2.375 @@ -1041,6 +1171,12 @@
   2.376              if (_excess[v] > 0 && _excess[v] <= delta)
   2.377                _active_nodes.push_back(v);
   2.378            }
   2.379 +
   2.380 +          // Global update heuristic
   2.381 +          if (relabel_cnt >= next_update_limit) {
   2.382 +            globalUpdate();
   2.383 +            next_update_limit += global_update_freq;
   2.384 +          }
   2.385          }
   2.386        }
   2.387      }
   2.388 @@ -1048,98 +1184,70 @@
   2.389      /// Execute the algorithm performing push and relabel operations
   2.390      void startPush() {
   2.391        // Paramters for heuristics
   2.392 -      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
   2.393 -      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
   2.394 +      const int EARLY_TERM_EPSILON_LIMIT = 1000;
   2.395 +      const double GLOBAL_UPDATE_FACTOR = 2.0;
   2.396  
   2.397 +      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
   2.398 +        (_res_node_num + _sup_node_num * _sup_node_num));
   2.399 +      int next_update_limit = global_update_freq;
   2.400 +
   2.401 +      int relabel_cnt = 0;
   2.402 +      
   2.403        // Perform cost scaling phases
   2.404        BoolVector hyper(_res_node_num, false);
   2.405 +      LargeCostVector hyper_cost(_res_node_num);
   2.406        for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
   2.407                                          1 : _epsilon / _alpha )
   2.408        {
   2.409 -        // "Early Termination" heuristic: use Bellman-Ford algorithm
   2.410 -        // to check if the current flow is optimal
   2.411 -        if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
   2.412 -          _arc_vec.clear();
   2.413 -          _cost_vec.clear();
   2.414 -          for (int j = 0; j != _res_arc_num; ++j) {
   2.415 -            if (_res_cap[j] > 0) {
   2.416 -              _arc_vec.push_back(IntPair(_source[j], _target[j]));
   2.417 -              _cost_vec.push_back(_cost[j] + 1);
   2.418 -            }
   2.419 -          }
   2.420 -          _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
   2.421 -
   2.422 -          BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
   2.423 -          bf.init(0);
   2.424 -          bool done = false;
   2.425 -          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
   2.426 -          for (int i = 0; i < K && !done; ++i)
   2.427 -            done = bf.processNextWeakRound();
   2.428 -          if (done) break;
   2.429 +        // Early termination heuristic
   2.430 +        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
   2.431 +          if (earlyTermination()) break;
   2.432          }
   2.433 -
   2.434 -        // Saturate arcs not satisfying the optimality condition
   2.435 -        for (int a = 0; a != _res_arc_num; ++a) {
   2.436 -          if (_res_cap[a] > 0 &&
   2.437 -              _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
   2.438 -            Value delta = _res_cap[a];
   2.439 -            _excess[_source[a]] -= delta;
   2.440 -            _excess[_target[a]] += delta;
   2.441 -            _res_cap[a] = 0;
   2.442 -            _res_cap[_reverse[a]] += delta;
   2.443 -          }
   2.444 -        }
   2.445 -
   2.446 -        // Find active nodes (i.e. nodes with positive excess)
   2.447 -        for (int u = 0; u != _res_node_num; ++u) {
   2.448 -          if (_excess[u] > 0) _active_nodes.push_back(u);
   2.449 -        }
   2.450 -
   2.451 -        // Initialize the next arcs
   2.452 -        for (int u = 0; u != _res_node_num; ++u) {
   2.453 -          _next_out[u] = _first_out[u];
   2.454 -        }
   2.455 +        
   2.456 +        // Initialize current phase
   2.457 +        initPhase();
   2.458  
   2.459          // Perform push and relabel operations
   2.460          while (_active_nodes.size() > 0) {
   2.461 -          LargeCost min_red_cost, rc;
   2.462 +          LargeCost min_red_cost, rc, pi_n;
   2.463            Value delta;
   2.464            int n, t, a, last_out = _res_arc_num;
   2.465  
   2.466 +        next_node:
   2.467            // Select an active node (FIFO selection)
   2.468 -        next_node:
   2.469            n = _active_nodes.front();
   2.470 -          last_out = _sum_supply < 0 ?
   2.471 -            _first_out[n+1] : _first_out[n+1] - 1;
   2.472 -
   2.473 +          last_out = _first_out[n+1];
   2.474 +          pi_n = _pi[n];
   2.475 +          
   2.476            // Perform push operations if there are admissible arcs
   2.477            if (_excess[n] > 0) {
   2.478              for (a = _next_out[n]; a != last_out; ++a) {
   2.479                if (_res_cap[a] > 0 &&
   2.480 -                  _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
   2.481 +                  _cost[a] + pi_n - _pi[_target[a]] < 0) {
   2.482                  delta = std::min(_res_cap[a], _excess[n]);
   2.483                  t = _target[a];
   2.484  
   2.485                  // Push-look-ahead heuristic
   2.486                  Value ahead = -_excess[t];
   2.487 -                int last_out_t = _sum_supply < 0 ?
   2.488 -                  _first_out[t+1] : _first_out[t+1] - 1;
   2.489 +                int last_out_t = _first_out[t+1];
   2.490 +                LargeCost pi_t = _pi[t];
   2.491                  for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
   2.492                    if (_res_cap[ta] > 0 && 
   2.493 -                      _cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0)
   2.494 +                      _cost[ta] + pi_t - _pi[_target[ta]] < 0)
   2.495                      ahead += _res_cap[ta];
   2.496                    if (ahead >= delta) break;
   2.497                  }
   2.498                  if (ahead < 0) ahead = 0;
   2.499  
   2.500                  // Push flow along the arc
   2.501 -                if (ahead < delta) {
   2.502 +                if (ahead < delta && !hyper[t]) {
   2.503                    _res_cap[a] -= ahead;
   2.504                    _res_cap[_reverse[a]] += ahead;
   2.505                    _excess[n] -= ahead;
   2.506                    _excess[t] += ahead;
   2.507                    _active_nodes.push_front(t);
   2.508                    hyper[t] = true;
   2.509 +                  hyper_cost[t] = _cost[a] + pi_n - pi_t;
   2.510                    _next_out[n] = a;
   2.511                    goto next_node;
   2.512                  } else {
   2.513 @@ -1162,18 +1270,18 @@
   2.514  
   2.515            // Relabel the node if it is still active (or hyper)
   2.516            if (_excess[n] > 0 || hyper[n]) {
   2.517 -            min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
   2.518 +             min_red_cost = hyper[n] ? -hyper_cost[n] :
   2.519 +               std::numeric_limits<LargeCost>::max();
   2.520              for (int a = _first_out[n]; a != last_out; ++a) {
   2.521 -              rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
   2.522 +              rc = _cost[a] + pi_n - _pi[_target[a]];
   2.523                if (_res_cap[a] > 0 && rc < min_red_cost) {
   2.524                  min_red_cost = rc;
   2.525                }
   2.526              }
   2.527              _pi[n] -= min_red_cost + _epsilon;
   2.528 +            _next_out[n] = _first_out[n];
   2.529              hyper[n] = false;
   2.530 -
   2.531 -            // Reset the next arc
   2.532 -            _next_out[n] = _first_out[n];
   2.533 +            ++relabel_cnt;
   2.534            }
   2.535          
   2.536            // Remove nodes that are not active nor hyper
   2.537 @@ -1183,6 +1291,14 @@
   2.538                    !hyper[_active_nodes.front()] ) {
   2.539              _active_nodes.pop_front();
   2.540            }
   2.541 +          
   2.542 +          // Global update heuristic
   2.543 +          if (relabel_cnt >= next_update_limit) {
   2.544 +            globalUpdate();
   2.545 +            for (int u = 0; u != _res_node_num; ++u)
   2.546 +              hyper[u] = false;
   2.547 +            next_update_limit += global_update_freq;
   2.548 +          }
   2.549          }
   2.550        }
   2.551      }
     3.1 --- a/lemon/cycle_canceling.h	Sun Feb 21 18:55:30 2010 +0100
     3.2 +++ b/lemon/cycle_canceling.h	Fri Feb 26 14:00:20 2010 +0100
     3.3 @@ -144,10 +144,11 @@
     3.4      TEMPLATE_DIGRAPH_TYPEDEFS(GR);
     3.5      
     3.6      typedef std::vector<int> IntVector;
     3.7 -    typedef std::vector<char> CharVector;
     3.8      typedef std::vector<double> DoubleVector;
     3.9      typedef std::vector<Value> ValueVector;
    3.10      typedef std::vector<Cost> CostVector;
    3.11 +    typedef std::vector<char> BoolVector;
    3.12 +    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
    3.13  
    3.14    private:
    3.15    
    3.16 @@ -198,7 +199,7 @@
    3.17      IntArcMap _arc_idf;
    3.18      IntArcMap _arc_idb;
    3.19      IntVector _first_out;
    3.20 -    CharVector _forward;
    3.21 +    BoolVector _forward;
    3.22      IntVector _source;
    3.23      IntVector _target;
    3.24      IntVector _reverse;
    3.25 @@ -962,8 +963,8 @@
    3.26        // Contruct auxiliary data vectors
    3.27        DoubleVector pi(_res_node_num, 0.0);
    3.28        IntVector level(_res_node_num);
    3.29 -      CharVector reached(_res_node_num);
    3.30 -      CharVector processed(_res_node_num);
    3.31 +      BoolVector reached(_res_node_num);
    3.32 +      BoolVector processed(_res_node_num);
    3.33        IntVector pred_node(_res_node_num);
    3.34        IntVector pred_arc(_res_node_num);
    3.35        std::vector<int> stack(_res_node_num);
     4.1 --- a/lemon/network_simplex.h	Sun Feb 21 18:55:30 2010 +0100
     4.2 +++ b/lemon/network_simplex.h	Fri Feb 26 14:00:20 2010 +0100
     4.3 @@ -164,9 +164,10 @@
     4.4      TEMPLATE_DIGRAPH_TYPEDEFS(GR);
     4.5  
     4.6      typedef std::vector<int> IntVector;
     4.7 -    typedef std::vector<char> CharVector;
     4.8      typedef std::vector<Value> ValueVector;
     4.9      typedef std::vector<Cost> CostVector;
    4.10 +    typedef std::vector<char> BoolVector;
    4.11 +    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
    4.12  
    4.13      // State constants for arcs
    4.14      enum ArcStateEnum {
    4.15 @@ -213,8 +214,8 @@
    4.16      IntVector _succ_num;
    4.17      IntVector _last_succ;
    4.18      IntVector _dirty_revs;
    4.19 -    CharVector _forward;
    4.20 -    CharVector _state;
    4.21 +    BoolVector _forward;
    4.22 +    BoolVector _state;
    4.23      int _root;
    4.24  
    4.25      // Temporary data used in the current pivot iteration
    4.26 @@ -245,7 +246,7 @@
    4.27        const IntVector  &_source;
    4.28        const IntVector  &_target;
    4.29        const CostVector &_cost;
    4.30 -      const CharVector &_state;
    4.31 +      const BoolVector &_state;
    4.32        const CostVector &_pi;
    4.33        int &_in_arc;
    4.34        int _search_arc_num;
    4.35 @@ -266,7 +267,7 @@
    4.36        // Find next entering arc
    4.37        bool findEnteringArc() {
    4.38          Cost c;
    4.39 -        for (int e = _next_arc; e < _search_arc_num; ++e) {
    4.40 +        for (int e = _next_arc; e != _search_arc_num; ++e) {
    4.41            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    4.42            if (c < 0) {
    4.43              _in_arc = e;
    4.44 @@ -274,7 +275,7 @@
    4.45              return true;
    4.46            }
    4.47          }
    4.48 -        for (int e = 0; e < _next_arc; ++e) {
    4.49 +        for (int e = 0; e != _next_arc; ++e) {
    4.50            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    4.51            if (c < 0) {
    4.52              _in_arc = e;
    4.53 @@ -297,7 +298,7 @@
    4.54        const IntVector  &_source;
    4.55        const IntVector  &_target;
    4.56        const CostVector &_cost;
    4.57 -      const CharVector &_state;
    4.58 +      const BoolVector &_state;
    4.59        const CostVector &_pi;
    4.60        int &_in_arc;
    4.61        int _search_arc_num;
    4.62 @@ -314,7 +315,7 @@
    4.63        // Find next entering arc
    4.64        bool findEnteringArc() {
    4.65          Cost c, min = 0;
    4.66 -        for (int e = 0; e < _search_arc_num; ++e) {
    4.67 +        for (int e = 0; e != _search_arc_num; ++e) {
    4.68            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    4.69            if (c < min) {
    4.70              min = c;
    4.71 @@ -336,7 +337,7 @@
    4.72        const IntVector  &_source;
    4.73        const IntVector  &_target;
    4.74        const CostVector &_cost;
    4.75 -      const CharVector &_state;
    4.76 +      const BoolVector &_state;
    4.77        const CostVector &_pi;
    4.78        int &_in_arc;
    4.79        int _search_arc_num;
    4.80 @@ -355,7 +356,7 @@
    4.81          _next_arc(0)
    4.82        {
    4.83          // The main parameters of the pivot rule
    4.84 -        const double BLOCK_SIZE_FACTOR = 0.5;
    4.85 +        const double BLOCK_SIZE_FACTOR = 1.0;
    4.86          const int MIN_BLOCK_SIZE = 10;
    4.87  
    4.88          _block_size = std::max( int(BLOCK_SIZE_FACTOR *
    4.89 @@ -368,7 +369,7 @@
    4.90          Cost c, min = 0;
    4.91          int cnt = _block_size;
    4.92          int e;
    4.93 -        for (e = _next_arc; e < _search_arc_num; ++e) {
    4.94 +        for (e = _next_arc; e != _search_arc_num; ++e) {
    4.95            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    4.96            if (c < min) {
    4.97              min = c;
    4.98 @@ -379,7 +380,7 @@
    4.99              cnt = _block_size;
   4.100            }
   4.101          }
   4.102 -        for (e = 0; e < _next_arc; ++e) {
   4.103 +        for (e = 0; e != _next_arc; ++e) {
   4.104            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   4.105            if (c < min) {
   4.106              min = c;
   4.107 @@ -409,7 +410,7 @@
   4.108        const IntVector  &_source;
   4.109        const IntVector  &_target;
   4.110        const CostVector &_cost;
   4.111 -      const CharVector &_state;
   4.112 +      const BoolVector &_state;
   4.113        const CostVector &_pi;
   4.114        int &_in_arc;
   4.115        int _search_arc_num;
   4.116 @@ -470,7 +471,7 @@
   4.117          // Major iteration: build a new candidate list
   4.118          min = 0;
   4.119          _curr_length = 0;
   4.120 -        for (e = _next_arc; e < _search_arc_num; ++e) {
   4.121 +        for (e = _next_arc; e != _search_arc_num; ++e) {
   4.122            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   4.123            if (c < 0) {
   4.124              _candidates[_curr_length++] = e;
   4.125 @@ -481,7 +482,7 @@
   4.126              if (_curr_length == _list_length) goto search_end;
   4.127            }
   4.128          }
   4.129 -        for (e = 0; e < _next_arc; ++e) {
   4.130 +        for (e = 0; e != _next_arc; ++e) {
   4.131            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   4.132            if (c < 0) {
   4.133              _candidates[_curr_length++] = e;
   4.134 @@ -512,7 +513,7 @@
   4.135        const IntVector  &_source;
   4.136        const IntVector  &_target;
   4.137        const CostVector &_cost;
   4.138 -      const CharVector &_state;
   4.139 +      const BoolVector &_state;
   4.140        const CostVector &_pi;
   4.141        int &_in_arc;
   4.142        int _search_arc_num;
   4.143 @@ -565,7 +566,7 @@
   4.144        bool findEnteringArc() {
   4.145          // Check the current candidate list
   4.146          int e;
   4.147 -        for (int i = 0; i < _curr_length; ++i) {
   4.148 +        for (int i = 0; i != _curr_length; ++i) {
   4.149            e = _candidates[i];
   4.150            _cand_cost[e] = _state[e] *
   4.151              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   4.152 @@ -578,7 +579,7 @@
   4.153          int cnt = _block_size;
   4.154          int limit = _head_length;
   4.155  
   4.156 -        for (e = _next_arc; e < _search_arc_num; ++e) {
   4.157 +        for (e = _next_arc; e != _search_arc_num; ++e) {
   4.158            _cand_cost[e] = _state[e] *
   4.159              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   4.160            if (_cand_cost[e] < 0) {
   4.161 @@ -590,7 +591,7 @@
   4.162              cnt = _block_size;
   4.163            }
   4.164          }
   4.165 -        for (e = 0; e < _next_arc; ++e) {
   4.166 +        for (e = 0; e != _next_arc; ++e) {
   4.167            _cand_cost[e] = _state[e] *
   4.168              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   4.169            if (_cand_cost[e] < 0) {
   4.170 @@ -1360,7 +1361,7 @@
   4.171        }
   4.172  
   4.173        // Update _rev_thread using the new _thread values
   4.174 -      for (int i = 0; i < int(_dirty_revs.size()); ++i) {
   4.175 +      for (int i = 0; i != int(_dirty_revs.size()); ++i) {
   4.176          u = _dirty_revs[i];
   4.177          _rev_thread[_thread[u]] = u;
   4.178        }
   4.179 @@ -1432,6 +1433,100 @@
   4.180        }
   4.181      }
   4.182  
   4.183 +    // Heuristic initial pivots
   4.184 +    bool initialPivots() {
   4.185 +      Value curr, total = 0;
   4.186 +      std::vector<Node> supply_nodes, demand_nodes;
   4.187 +      for (NodeIt u(_graph); u != INVALID; ++u) {
   4.188 +        curr = _supply[_node_id[u]];
   4.189 +        if (curr > 0) {
   4.190 +          total += curr;
   4.191 +          supply_nodes.push_back(u);
   4.192 +        }
   4.193 +        else if (curr < 0) {
   4.194 +          demand_nodes.push_back(u);
   4.195 +        }
   4.196 +      }
   4.197 +      if (_sum_supply > 0) total -= _sum_supply;
   4.198 +      if (total <= 0) return true;
   4.199 +
   4.200 +      IntVector arc_vector;
   4.201 +      if (_sum_supply >= 0) {
   4.202 +        if (supply_nodes.size() == 1 && demand_nodes.size() == 1) {
   4.203 +          // Perform a reverse graph search from the sink to the source
   4.204 +          typename GR::template NodeMap<bool> reached(_graph, false);
   4.205 +          Node s = supply_nodes[0], t = demand_nodes[0];
   4.206 +          std::vector<Node> stack;
   4.207 +          reached[t] = true;
   4.208 +          stack.push_back(t);
   4.209 +          while (!stack.empty()) {
   4.210 +            Node u, v = stack.back();
   4.211 +            stack.pop_back();
   4.212 +            if (v == s) break;
   4.213 +            for (InArcIt a(_graph, v); a != INVALID; ++a) {
   4.214 +              if (reached[u = _graph.source(a)]) continue;
   4.215 +              int j = _arc_id[a];
   4.216 +              if (_cap[j] >= total) {
   4.217 +                arc_vector.push_back(j);
   4.218 +                reached[u] = true;
   4.219 +                stack.push_back(u);
   4.220 +              }
   4.221 +            }
   4.222 +          }
   4.223 +        } else {
   4.224 +          // Find the min. cost incomming arc for each demand node
   4.225 +          for (int i = 0; i != int(demand_nodes.size()); ++i) {
   4.226 +            Node v = demand_nodes[i];
   4.227 +            Cost c, min_cost = std::numeric_limits<Cost>::max();
   4.228 +            Arc min_arc = INVALID;
   4.229 +            for (InArcIt a(_graph, v); a != INVALID; ++a) {
   4.230 +              c = _cost[_arc_id[a]];
   4.231 +              if (c < min_cost) {
   4.232 +                min_cost = c;
   4.233 +                min_arc = a;
   4.234 +              }
   4.235 +            }
   4.236 +            if (min_arc != INVALID) {
   4.237 +              arc_vector.push_back(_arc_id[min_arc]);
   4.238 +            }
   4.239 +          }
   4.240 +        }
   4.241 +      } else {
   4.242 +        // Find the min. cost outgoing arc for each supply node
   4.243 +        for (int i = 0; i != int(supply_nodes.size()); ++i) {
   4.244 +          Node u = supply_nodes[i];
   4.245 +          Cost c, min_cost = std::numeric_limits<Cost>::max();
   4.246 +          Arc min_arc = INVALID;
   4.247 +          for (OutArcIt a(_graph, u); a != INVALID; ++a) {
   4.248 +            c = _cost[_arc_id[a]];
   4.249 +            if (c < min_cost) {
   4.250 +              min_cost = c;
   4.251 +              min_arc = a;
   4.252 +            }
   4.253 +          }
   4.254 +          if (min_arc != INVALID) {
   4.255 +            arc_vector.push_back(_arc_id[min_arc]);
   4.256 +          }
   4.257 +        }
   4.258 +      }
   4.259 +
   4.260 +      // Perform heuristic initial pivots
   4.261 +      for (int i = 0; i != int(arc_vector.size()); ++i) {
   4.262 +        in_arc = arc_vector[i];
   4.263 +        if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] -
   4.264 +            _pi[_target[in_arc]]) >= 0) continue;
   4.265 +        findJoinNode();
   4.266 +        bool change = findLeavingArc();
   4.267 +        if (delta >= MAX) return false;
   4.268 +        changeFlow(change);
   4.269 +        if (change) {
   4.270 +          updateTreeStructure();
   4.271 +          updatePotential();
   4.272 +        }
   4.273 +      }
   4.274 +      return true;
   4.275 +    }
   4.276 +
   4.277      // Execute the algorithm
   4.278      ProblemType start(PivotRule pivot_rule) {
   4.279        // Select the pivot rule implementation
   4.280 @@ -1454,6 +1549,9 @@
   4.281      ProblemType start() {
   4.282        PivotRuleImpl pivot(*this);
   4.283  
   4.284 +      // Perform heuristic initial pivots
   4.285 +      if (!initialPivots()) return UNBOUNDED;
   4.286 +
   4.287        // Execute the Network Simplex algorithm
   4.288        while (pivot.findEnteringArc()) {
   4.289          findJoinNode();