lemon/cost_scaling.h
changeset 839 f3bc4e9b5f3a
parent 821 072ec8120958
child 840 2914b6f0fde0
     1.1 --- a/lemon/cost_scaling.h	Tue Feb 09 23:29:51 2010 +0100
     1.2 +++ b/lemon/cost_scaling.h	Sat Feb 20 18:39:03 2010 +0100
     1.3 @@ -197,10 +197,11 @@
     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<LargeCost> LargeCostVector;
    1.11 +    typedef std::vector<char> BoolVector;
    1.12 +    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
    1.13  
    1.14    private:
    1.15    
    1.16 @@ -244,6 +245,7 @@
    1.17      // Parameters of the problem
    1.18      bool _have_lower;
    1.19      Value _sum_supply;
    1.20 +    int _sup_node_num;
    1.21  
    1.22      // Data structures for storing the digraph
    1.23      IntNodeMap _node_id;
    1.24 @@ -272,6 +274,12 @@
    1.25      LargeCost _epsilon;
    1.26      int _alpha;
    1.27  
    1.28 +    IntVector _buckets;
    1.29 +    IntVector _bucket_next;
    1.30 +    IntVector _bucket_prev;
    1.31 +    IntVector _rank;
    1.32 +    int _max_rank;
    1.33 +  
    1.34      // Data for a StaticDigraph structure
    1.35      typedef std::pair<int, int> IntPair;
    1.36      StaticDigraph _sgr;
    1.37 @@ -802,6 +810,11 @@
    1.38          }
    1.39        }
    1.40  
    1.41 +      _sup_node_num = 0;
    1.42 +      for (NodeIt n(_graph); n != INVALID; ++n) {
    1.43 +        if (sup[n] > 0) ++_sup_node_num;
    1.44 +      }
    1.45 +
    1.46        // Find a feasible flow using Circulation
    1.47        Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
    1.48          circ(_graph, low, cap, sup);
    1.49 @@ -836,7 +849,7 @@
    1.50          }
    1.51          for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
    1.52            int ra = _reverse[a];
    1.53 -          _res_cap[a] = 1;
    1.54 +          _res_cap[a] = 0;
    1.55            _res_cap[ra] = 0;
    1.56            _cost[a] = 0;
    1.57            _cost[ra] = 0;
    1.58 @@ -850,7 +863,14 @@
    1.59      void start(Method method) {
    1.60        // Maximum path length for partial augment
    1.61        const int MAX_PATH_LENGTH = 4;
    1.62 -      
    1.63 +
    1.64 +      // Initialize data structures for buckets      
    1.65 +      _max_rank = _alpha * _res_node_num;
    1.66 +      _buckets.resize(_max_rank);
    1.67 +      _bucket_next.resize(_res_node_num + 1);
    1.68 +      _bucket_prev.resize(_res_node_num + 1);
    1.69 +      _rank.resize(_res_node_num + 1);
    1.70 +  
    1.71        // Execute the algorithm
    1.72        switch (method) {
    1.73          case PUSH:
    1.74 @@ -889,63 +909,175 @@
    1.75          }
    1.76        }
    1.77      }
    1.78 +    
    1.79 +    // Initialize a cost scaling phase
    1.80 +    void initPhase() {
    1.81 +      // Saturate arcs not satisfying the optimality condition
    1.82 +      for (int u = 0; u != _res_node_num; ++u) {
    1.83 +        int last_out = _first_out[u+1];
    1.84 +        LargeCost pi_u = _pi[u];
    1.85 +        for (int a = _first_out[u]; a != last_out; ++a) {
    1.86 +          int v = _target[a];
    1.87 +          if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) {
    1.88 +            Value delta = _res_cap[a];
    1.89 +            _excess[u] -= delta;
    1.90 +            _excess[v] += delta;
    1.91 +            _res_cap[a] = 0;
    1.92 +            _res_cap[_reverse[a]] += delta;
    1.93 +          }
    1.94 +        }
    1.95 +      }
    1.96 +      
    1.97 +      // Find active nodes (i.e. nodes with positive excess)
    1.98 +      for (int u = 0; u != _res_node_num; ++u) {
    1.99 +        if (_excess[u] > 0) _active_nodes.push_back(u);
   1.100 +      }
   1.101 +
   1.102 +      // Initialize the next arcs
   1.103 +      for (int u = 0; u != _res_node_num; ++u) {
   1.104 +        _next_out[u] = _first_out[u];
   1.105 +      }
   1.106 +    }
   1.107 +    
   1.108 +    // Early termination heuristic
   1.109 +    bool earlyTermination() {
   1.110 +      const double EARLY_TERM_FACTOR = 3.0;
   1.111 +
   1.112 +      // Build a static residual graph
   1.113 +      _arc_vec.clear();
   1.114 +      _cost_vec.clear();
   1.115 +      for (int j = 0; j != _res_arc_num; ++j) {
   1.116 +        if (_res_cap[j] > 0) {
   1.117 +          _arc_vec.push_back(IntPair(_source[j], _target[j]));
   1.118 +          _cost_vec.push_back(_cost[j] + 1);
   1.119 +        }
   1.120 +      }
   1.121 +      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
   1.122 +
   1.123 +      // Run Bellman-Ford algorithm to check if the current flow is optimal
   1.124 +      BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
   1.125 +      bf.init(0);
   1.126 +      bool done = false;
   1.127 +      int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num)));
   1.128 +      for (int i = 0; i < K && !done; ++i) {
   1.129 +        done = bf.processNextWeakRound();
   1.130 +      }
   1.131 +      return done;
   1.132 +    }
   1.133 +
   1.134 +    // Global potential update heuristic
   1.135 +    void globalUpdate() {
   1.136 +      int bucket_end = _root + 1;
   1.137 +    
   1.138 +      // Initialize buckets
   1.139 +      for (int r = 0; r != _max_rank; ++r) {
   1.140 +        _buckets[r] = bucket_end;
   1.141 +      }
   1.142 +      Value total_excess = 0;
   1.143 +      for (int i = 0; i != _res_node_num; ++i) {
   1.144 +        if (_excess[i] < 0) {
   1.145 +          _rank[i] = 0;
   1.146 +          _bucket_next[i] = _buckets[0];
   1.147 +          _bucket_prev[_buckets[0]] = i;
   1.148 +          _buckets[0] = i;
   1.149 +        } else {
   1.150 +          total_excess += _excess[i];
   1.151 +          _rank[i] = _max_rank;
   1.152 +        }
   1.153 +      }
   1.154 +      if (total_excess == 0) return;
   1.155 +
   1.156 +      // Search the buckets
   1.157 +      int r = 0;
   1.158 +      for ( ; r != _max_rank; ++r) {
   1.159 +        while (_buckets[r] != bucket_end) {
   1.160 +          // Remove the first node from the current bucket
   1.161 +          int u = _buckets[r];
   1.162 +          _buckets[r] = _bucket_next[u];
   1.163 +          
   1.164 +          // Search the incomming arcs of u
   1.165 +          LargeCost pi_u = _pi[u];
   1.166 +          int last_out = _first_out[u+1];
   1.167 +          for (int a = _first_out[u]; a != last_out; ++a) {
   1.168 +            int ra = _reverse[a];
   1.169 +            if (_res_cap[ra] > 0) {
   1.170 +              int v = _source[ra];
   1.171 +              int old_rank_v = _rank[v];
   1.172 +              if (r < old_rank_v) {
   1.173 +                // Compute the new rank of v
   1.174 +                LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
   1.175 +                int new_rank_v = old_rank_v;
   1.176 +                if (nrc < LargeCost(_max_rank))
   1.177 +                  new_rank_v = r + 1 + int(nrc);
   1.178 +                  
   1.179 +                // Change the rank of v
   1.180 +                if (new_rank_v < old_rank_v) {
   1.181 +                  _rank[v] = new_rank_v;
   1.182 +                  _next_out[v] = _first_out[v];
   1.183 +                  
   1.184 +                  // Remove v from its old bucket
   1.185 +                  if (old_rank_v < _max_rank) {
   1.186 +                    if (_buckets[old_rank_v] == v) {
   1.187 +                      _buckets[old_rank_v] = _bucket_next[v];
   1.188 +                    } else {
   1.189 +                      _bucket_next[_bucket_prev[v]] = _bucket_next[v];
   1.190 +                      _bucket_prev[_bucket_next[v]] = _bucket_prev[v];
   1.191 +                    }
   1.192 +                  }
   1.193 +                  
   1.194 +                  // Insert v to its new bucket
   1.195 +                  _bucket_next[v] = _buckets[new_rank_v];
   1.196 +                  _bucket_prev[_buckets[new_rank_v]] = v;
   1.197 +                  _buckets[new_rank_v] = v;
   1.198 +                }
   1.199 +              }
   1.200 +            }
   1.201 +          }
   1.202 +
   1.203 +          // Finish search if there are no more active nodes
   1.204 +          if (_excess[u] > 0) {
   1.205 +            total_excess -= _excess[u];
   1.206 +            if (total_excess <= 0) break;
   1.207 +          }
   1.208 +        }
   1.209 +        if (total_excess <= 0) break;
   1.210 +      }
   1.211 +      
   1.212 +      // Relabel nodes
   1.213 +      for (int u = 0; u != _res_node_num; ++u) {
   1.214 +        int k = std::min(_rank[u], r);
   1.215 +        if (k > 0) {
   1.216 +          _pi[u] -= _epsilon * k;
   1.217 +          _next_out[u] = _first_out[u];
   1.218 +        }
   1.219 +      }
   1.220 +    }
   1.221  
   1.222      /// Execute the algorithm performing augment and relabel operations
   1.223      void startAugment(int max_length = std::numeric_limits<int>::max()) {
   1.224        // Paramters for heuristics
   1.225 -      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
   1.226 -      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
   1.227 +      const int EARLY_TERM_EPSILON_LIMIT = 1000;
   1.228 +      const double GLOBAL_UPDATE_FACTOR = 3.0;
   1.229  
   1.230 +      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
   1.231 +        (_res_node_num + _sup_node_num * _sup_node_num));
   1.232 +      int next_update_limit = global_update_freq;
   1.233 +      
   1.234 +      int relabel_cnt = 0;
   1.235 +      
   1.236        // Perform cost scaling phases
   1.237 -      IntVector pred_arc(_res_node_num);
   1.238 -      std::vector<int> path_nodes;
   1.239 +      std::vector<int> path;
   1.240        for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
   1.241                                          1 : _epsilon / _alpha )
   1.242        {
   1.243 -        // "Early Termination" heuristic: use Bellman-Ford algorithm
   1.244 -        // to check if the current flow is optimal
   1.245 -        if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
   1.246 -          _arc_vec.clear();
   1.247 -          _cost_vec.clear();
   1.248 -          for (int j = 0; j != _res_arc_num; ++j) {
   1.249 -            if (_res_cap[j] > 0) {
   1.250 -              _arc_vec.push_back(IntPair(_source[j], _target[j]));
   1.251 -              _cost_vec.push_back(_cost[j] + 1);
   1.252 -            }
   1.253 -          }
   1.254 -          _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
   1.255 -
   1.256 -          BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
   1.257 -          bf.init(0);
   1.258 -          bool done = false;
   1.259 -          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
   1.260 -          for (int i = 0; i < K && !done; ++i)
   1.261 -            done = bf.processNextWeakRound();
   1.262 -          if (done) break;
   1.263 -        }
   1.264 -
   1.265 -        // Saturate arcs not satisfying the optimality condition
   1.266 -        for (int a = 0; a != _res_arc_num; ++a) {
   1.267 -          if (_res_cap[a] > 0 &&
   1.268 -              _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
   1.269 -            Value delta = _res_cap[a];
   1.270 -            _excess[_source[a]] -= delta;
   1.271 -            _excess[_target[a]] += delta;
   1.272 -            _res_cap[a] = 0;
   1.273 -            _res_cap[_reverse[a]] += delta;
   1.274 -          }
   1.275 +        // Early termination heuristic
   1.276 +        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
   1.277 +          if (earlyTermination()) break;
   1.278          }
   1.279          
   1.280 -        // Find active nodes (i.e. nodes with positive excess)
   1.281 -        for (int u = 0; u != _res_node_num; ++u) {
   1.282 -          if (_excess[u] > 0) _active_nodes.push_back(u);
   1.283 -        }
   1.284 -
   1.285 -        // Initialize the next arcs
   1.286 -        for (int u = 0; u != _res_node_num; ++u) {
   1.287 -          _next_out[u] = _first_out[u];
   1.288 -        }
   1.289 -
   1.290 +        // Initialize current phase
   1.291 +        initPhase();
   1.292 +        
   1.293          // Perform partial augment and relabel operations
   1.294          while (true) {
   1.295            // Select an active node (FIFO selection)
   1.296 @@ -955,46 +1087,44 @@
   1.297            }
   1.298            if (_active_nodes.size() == 0) break;
   1.299            int start = _active_nodes.front();
   1.300 -          path_nodes.clear();
   1.301 -          path_nodes.push_back(start);
   1.302  
   1.303            // Find an augmenting path from the start node
   1.304 +          path.clear();
   1.305            int tip = start;
   1.306 -          while (_excess[tip] >= 0 &&
   1.307 -                 int(path_nodes.size()) <= max_length) {
   1.308 +          while (_excess[tip] >= 0 && int(path.size()) < max_length) {
   1.309              int u;
   1.310 -            LargeCost min_red_cost, rc;
   1.311 -            int last_out = _sum_supply < 0 ?
   1.312 -              _first_out[tip+1] : _first_out[tip+1] - 1;
   1.313 +            LargeCost min_red_cost, rc, pi_tip = _pi[tip];
   1.314 +            int last_out = _first_out[tip+1];
   1.315              for (int a = _next_out[tip]; a != last_out; ++a) {
   1.316 -              if (_res_cap[a] > 0 &&
   1.317 -                  _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
   1.318 -                u = _target[a];
   1.319 -                pred_arc[u] = a;
   1.320 +              u = _target[a];
   1.321 +              if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) {
   1.322 +                path.push_back(a);
   1.323                  _next_out[tip] = a;
   1.324                  tip = u;
   1.325 -                path_nodes.push_back(tip);
   1.326                  goto next_step;
   1.327                }
   1.328              }
   1.329  
   1.330              // Relabel tip node
   1.331 -            min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
   1.332 +            min_red_cost = std::numeric_limits<LargeCost>::max();
   1.333 +            if (tip != start) {
   1.334 +              int ra = _reverse[path.back()];
   1.335 +              min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]];
   1.336 +            }
   1.337              for (int a = _first_out[tip]; a != last_out; ++a) {
   1.338 -              rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
   1.339 +              rc = _cost[a] + pi_tip - _pi[_target[a]];
   1.340                if (_res_cap[a] > 0 && rc < min_red_cost) {
   1.341                  min_red_cost = rc;
   1.342                }
   1.343              }
   1.344              _pi[tip] -= min_red_cost + _epsilon;
   1.345 -
   1.346 -            // Reset the next arc of tip
   1.347              _next_out[tip] = _first_out[tip];
   1.348 +            ++relabel_cnt;
   1.349  
   1.350              // Step back
   1.351              if (tip != start) {
   1.352 -              path_nodes.pop_back();
   1.353 -              tip = path_nodes.back();
   1.354 +              tip = _source[path.back()];
   1.355 +              path.pop_back();
   1.356              }
   1.357  
   1.358            next_step: ;
   1.359 @@ -1002,11 +1132,11 @@
   1.360  
   1.361            // Augment along the found path (as much flow as possible)
   1.362            Value delta;
   1.363 -          int u, v = path_nodes.front(), pa;
   1.364 -          for (int i = 1; i < int(path_nodes.size()); ++i) {
   1.365 +          int pa, u, v = start;
   1.366 +          for (int i = 0; i != int(path.size()); ++i) {
   1.367 +            pa = path[i];
   1.368              u = v;
   1.369 -            v = path_nodes[i];
   1.370 -            pa = pred_arc[v];
   1.371 +            v = _target[pa];
   1.372              delta = std::min(_res_cap[pa], _excess[u]);
   1.373              _res_cap[pa] -= delta;
   1.374              _res_cap[_reverse[pa]] += delta;
   1.375 @@ -1015,6 +1145,12 @@
   1.376              if (_excess[v] > 0 && _excess[v] <= delta)
   1.377                _active_nodes.push_back(v);
   1.378            }
   1.379 +
   1.380 +          // Global update heuristic
   1.381 +          if (relabel_cnt >= next_update_limit) {
   1.382 +            globalUpdate();
   1.383 +            next_update_limit += global_update_freq;
   1.384 +          }
   1.385          }
   1.386        }
   1.387      }
   1.388 @@ -1022,98 +1158,70 @@
   1.389      /// Execute the algorithm performing push and relabel operations
   1.390      void startPush() {
   1.391        // Paramters for heuristics
   1.392 -      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
   1.393 -      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
   1.394 +      const int EARLY_TERM_EPSILON_LIMIT = 1000;
   1.395 +      const double GLOBAL_UPDATE_FACTOR = 2.0;
   1.396  
   1.397 +      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
   1.398 +        (_res_node_num + _sup_node_num * _sup_node_num));
   1.399 +      int next_update_limit = global_update_freq;
   1.400 +
   1.401 +      int relabel_cnt = 0;
   1.402 +      
   1.403        // Perform cost scaling phases
   1.404        BoolVector hyper(_res_node_num, false);
   1.405 +      LargeCostVector hyper_cost(_res_node_num);
   1.406        for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
   1.407                                          1 : _epsilon / _alpha )
   1.408        {
   1.409 -        // "Early Termination" heuristic: use Bellman-Ford algorithm
   1.410 -        // to check if the current flow is optimal
   1.411 -        if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
   1.412 -          _arc_vec.clear();
   1.413 -          _cost_vec.clear();
   1.414 -          for (int j = 0; j != _res_arc_num; ++j) {
   1.415 -            if (_res_cap[j] > 0) {
   1.416 -              _arc_vec.push_back(IntPair(_source[j], _target[j]));
   1.417 -              _cost_vec.push_back(_cost[j] + 1);
   1.418 -            }
   1.419 -          }
   1.420 -          _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
   1.421 -
   1.422 -          BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
   1.423 -          bf.init(0);
   1.424 -          bool done = false;
   1.425 -          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
   1.426 -          for (int i = 0; i < K && !done; ++i)
   1.427 -            done = bf.processNextWeakRound();
   1.428 -          if (done) break;
   1.429 +        // Early termination heuristic
   1.430 +        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
   1.431 +          if (earlyTermination()) break;
   1.432          }
   1.433 -
   1.434 -        // Saturate arcs not satisfying the optimality condition
   1.435 -        for (int a = 0; a != _res_arc_num; ++a) {
   1.436 -          if (_res_cap[a] > 0 &&
   1.437 -              _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
   1.438 -            Value delta = _res_cap[a];
   1.439 -            _excess[_source[a]] -= delta;
   1.440 -            _excess[_target[a]] += delta;
   1.441 -            _res_cap[a] = 0;
   1.442 -            _res_cap[_reverse[a]] += delta;
   1.443 -          }
   1.444 -        }
   1.445 -
   1.446 -        // Find active nodes (i.e. nodes with positive excess)
   1.447 -        for (int u = 0; u != _res_node_num; ++u) {
   1.448 -          if (_excess[u] > 0) _active_nodes.push_back(u);
   1.449 -        }
   1.450 -
   1.451 -        // Initialize the next arcs
   1.452 -        for (int u = 0; u != _res_node_num; ++u) {
   1.453 -          _next_out[u] = _first_out[u];
   1.454 -        }
   1.455 +        
   1.456 +        // Initialize current phase
   1.457 +        initPhase();
   1.458  
   1.459          // Perform push and relabel operations
   1.460          while (_active_nodes.size() > 0) {
   1.461 -          LargeCost min_red_cost, rc;
   1.462 +          LargeCost min_red_cost, rc, pi_n;
   1.463            Value delta;
   1.464            int n, t, a, last_out = _res_arc_num;
   1.465  
   1.466 +        next_node:
   1.467            // Select an active node (FIFO selection)
   1.468 -        next_node:
   1.469            n = _active_nodes.front();
   1.470 -          last_out = _sum_supply < 0 ?
   1.471 -            _first_out[n+1] : _first_out[n+1] - 1;
   1.472 -
   1.473 +          last_out = _first_out[n+1];
   1.474 +          pi_n = _pi[n];
   1.475 +          
   1.476            // Perform push operations if there are admissible arcs
   1.477            if (_excess[n] > 0) {
   1.478              for (a = _next_out[n]; a != last_out; ++a) {
   1.479                if (_res_cap[a] > 0 &&
   1.480 -                  _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
   1.481 +                  _cost[a] + pi_n - _pi[_target[a]] < 0) {
   1.482                  delta = std::min(_res_cap[a], _excess[n]);
   1.483                  t = _target[a];
   1.484  
   1.485                  // Push-look-ahead heuristic
   1.486                  Value ahead = -_excess[t];
   1.487 -                int last_out_t = _sum_supply < 0 ?
   1.488 -                  _first_out[t+1] : _first_out[t+1] - 1;
   1.489 +                int last_out_t = _first_out[t+1];
   1.490 +                LargeCost pi_t = _pi[t];
   1.491                  for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
   1.492                    if (_res_cap[ta] > 0 && 
   1.493 -                      _cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0)
   1.494 +                      _cost[ta] + pi_t - _pi[_target[ta]] < 0)
   1.495                      ahead += _res_cap[ta];
   1.496                    if (ahead >= delta) break;
   1.497                  }
   1.498                  if (ahead < 0) ahead = 0;
   1.499  
   1.500                  // Push flow along the arc
   1.501 -                if (ahead < delta) {
   1.502 +                if (ahead < delta && !hyper[t]) {
   1.503                    _res_cap[a] -= ahead;
   1.504                    _res_cap[_reverse[a]] += ahead;
   1.505                    _excess[n] -= ahead;
   1.506                    _excess[t] += ahead;
   1.507                    _active_nodes.push_front(t);
   1.508                    hyper[t] = true;
   1.509 +                  hyper_cost[t] = _cost[a] + pi_n - pi_t;
   1.510                    _next_out[n] = a;
   1.511                    goto next_node;
   1.512                  } else {
   1.513 @@ -1136,18 +1244,18 @@
   1.514  
   1.515            // Relabel the node if it is still active (or hyper)
   1.516            if (_excess[n] > 0 || hyper[n]) {
   1.517 -            min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
   1.518 +             min_red_cost = hyper[n] ? -hyper_cost[n] :
   1.519 +               std::numeric_limits<LargeCost>::max();
   1.520              for (int a = _first_out[n]; a != last_out; ++a) {
   1.521 -              rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
   1.522 +              rc = _cost[a] + pi_n - _pi[_target[a]];
   1.523                if (_res_cap[a] > 0 && rc < min_red_cost) {
   1.524                  min_red_cost = rc;
   1.525                }
   1.526              }
   1.527              _pi[n] -= min_red_cost + _epsilon;
   1.528 +            _next_out[n] = _first_out[n];
   1.529              hyper[n] = false;
   1.530 -
   1.531 -            // Reset the next arc
   1.532 -            _next_out[n] = _first_out[n];
   1.533 +            ++relabel_cnt;
   1.534            }
   1.535          
   1.536            // Remove nodes that are not active nor hyper
   1.537 @@ -1157,6 +1265,14 @@
   1.538                    !hyper[_active_nodes.front()] ) {
   1.539              _active_nodes.pop_front();
   1.540            }
   1.541 +          
   1.542 +          // Global update heuristic
   1.543 +          if (relabel_cnt >= next_update_limit) {
   1.544 +            globalUpdate();
   1.545 +            for (int u = 0; u != _res_node_num; ++u)
   1.546 +              hyper[u] = false;
   1.547 +            next_update_limit += global_update_freq;
   1.548 +          }
   1.549          }
   1.550        }
   1.551      }