Fix and improve refine methods in CostScaling (#417)
authorPeter Kovacs <kpeter@inf.elte.hu>
Tue, 15 Mar 2011 19:16:20 +0100 (2011-03-15)
changeset 9356ea176638264
parent 934 fe283caf6414
child 936 ddd3c0d3d9bf
Fix and improve refine methods in CostScaling (#417)
lemon/cost_scaling.h
     1.1 --- a/lemon/cost_scaling.h	Tue Mar 15 17:59:57 2011 +0100
     1.2 +++ b/lemon/cost_scaling.h	Tue Mar 15 19:16:20 2011 +0100
     1.3 @@ -1103,16 +1103,15 @@
     1.4      void startAugment(int max_length) {
     1.5        // Paramters for heuristics
     1.6        const int EARLY_TERM_EPSILON_LIMIT = 1000;
     1.7 -      const double GLOBAL_UPDATE_FACTOR = 3.0;
     1.8 -
     1.9 -      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
    1.10 +      const double GLOBAL_UPDATE_FACTOR = 1.0;
    1.11 +      const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
    1.12          (_res_node_num + _sup_node_num * _sup_node_num));
    1.13 -      int next_update_limit = global_update_freq;
    1.14 -
    1.15 -      int relabel_cnt = 0;
    1.16 +      int next_global_update_limit = global_update_skip;
    1.17  
    1.18        // Perform cost scaling phases
    1.19 -      std::vector<int> path;
    1.20 +      IntVector path;
    1.21 +      BoolVector path_arc(_res_arc_num, false);
    1.22 +      int relabel_cnt = 0;
    1.23        for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
    1.24                                          1 : _epsilon / _alpha )
    1.25        {
    1.26 @@ -1135,32 +1134,45 @@
    1.27            int start = _active_nodes.front();
    1.28  
    1.29            // Find an augmenting path from the start node
    1.30 -          path.clear();
    1.31            int tip = start;
    1.32 -          while (_excess[tip] >= 0 && int(path.size()) < max_length) {
    1.33 +          while (int(path.size()) < max_length && _excess[tip] >= 0) {
    1.34              int u;
    1.35 -            LargeCost min_red_cost, rc, pi_tip = _pi[tip];
    1.36 +            LargeCost rc, min_red_cost = std::numeric_limits<LargeCost>::max();
    1.37 +            LargeCost pi_tip = _pi[tip];
    1.38              int last_out = _first_out[tip+1];
    1.39              for (int a = _next_out[tip]; a != last_out; ++a) {
    1.40 -              u = _target[a];
    1.41 -              if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) {
    1.42 -                path.push_back(a);
    1.43 -                _next_out[tip] = a;
    1.44 -                tip = u;
    1.45 -                goto next_step;
    1.46 +              if (_res_cap[a] > 0) {
    1.47 +                u = _target[a];
    1.48 +                rc = _cost[a] + pi_tip - _pi[u];
    1.49 +                if (rc < 0) {
    1.50 +                  path.push_back(a);
    1.51 +                  _next_out[tip] = a;
    1.52 +                  if (path_arc[a]) {
    1.53 +                    goto augment;   // a cycle is found, stop path search
    1.54 +                  }
    1.55 +                  tip = u;
    1.56 +                  path_arc[a] = true;
    1.57 +                  goto next_step;
    1.58 +                }
    1.59 +                else if (rc < min_red_cost) {
    1.60 +                  min_red_cost = rc;
    1.61 +                }
    1.62                }
    1.63              }
    1.64  
    1.65              // Relabel tip node
    1.66 -            min_red_cost = std::numeric_limits<LargeCost>::max();
    1.67              if (tip != start) {
    1.68                int ra = _reverse[path.back()];
    1.69 -              min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]];
    1.70 +              min_red_cost =
    1.71 +                std::min(min_red_cost, _cost[ra] + pi_tip - _pi[_target[ra]]);
    1.72              }
    1.73 +            last_out = _next_out[tip];
    1.74              for (int a = _first_out[tip]; a != last_out; ++a) {
    1.75 -              rc = _cost[a] + pi_tip - _pi[_target[a]];
    1.76 -              if (_res_cap[a] > 0 && rc < min_red_cost) {
    1.77 -                min_red_cost = rc;
    1.78 +              if (_res_cap[a] > 0) {
    1.79 +                rc = _cost[a] + pi_tip - _pi[_target[a]];
    1.80 +                if (rc < min_red_cost) {
    1.81 +                  min_red_cost = rc;
    1.82 +                }
    1.83                }
    1.84              }
    1.85              _pi[tip] -= min_red_cost + _epsilon;
    1.86 @@ -1169,7 +1181,9 @@
    1.87  
    1.88              // Step back
    1.89              if (tip != start) {
    1.90 -              tip = _source[path.back()];
    1.91 +              int pa = path.back();
    1.92 +              path_arc[pa] = false;
    1.93 +              tip = _source[pa];
    1.94                path.pop_back();
    1.95              }
    1.96  
    1.97 @@ -1177,28 +1191,34 @@
    1.98            }
    1.99  
   1.100            // Augment along the found path (as much flow as possible)
   1.101 +        augment:
   1.102            Value delta;
   1.103            int pa, u, v = start;
   1.104            for (int i = 0; i != int(path.size()); ++i) {
   1.105              pa = path[i];
   1.106              u = v;
   1.107              v = _target[pa];
   1.108 +            path_arc[pa] = false;
   1.109              delta = std::min(_res_cap[pa], _excess[u]);
   1.110              _res_cap[pa] -= delta;
   1.111              _res_cap[_reverse[pa]] += delta;
   1.112              _excess[u] -= delta;
   1.113              _excess[v] += delta;
   1.114 -            if (_excess[v] > 0 && _excess[v] <= delta)
   1.115 +            if (_excess[v] > 0 && _excess[v] <= delta) {
   1.116                _active_nodes.push_back(v);
   1.117 +            }
   1.118            }
   1.119 +          path.clear();
   1.120  
   1.121            // Global update heuristic
   1.122 -          if (relabel_cnt >= next_update_limit) {
   1.123 +          if (relabel_cnt >= next_global_update_limit) {
   1.124              globalUpdate();
   1.125 -            next_update_limit += global_update_freq;
   1.126 +            next_global_update_limit += global_update_skip;
   1.127            }
   1.128          }
   1.129 +
   1.130        }
   1.131 +
   1.132      }
   1.133  
   1.134      /// Execute the algorithm performing push and relabel operations
   1.135 @@ -1207,15 +1227,14 @@
   1.136        const int EARLY_TERM_EPSILON_LIMIT = 1000;
   1.137        const double GLOBAL_UPDATE_FACTOR = 2.0;
   1.138  
   1.139 -      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
   1.140 +      const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
   1.141          (_res_node_num + _sup_node_num * _sup_node_num));
   1.142 -      int next_update_limit = global_update_freq;
   1.143 -
   1.144 -      int relabel_cnt = 0;
   1.145 +      int next_global_update_limit = global_update_skip;
   1.146  
   1.147        // Perform cost scaling phases
   1.148        BoolVector hyper(_res_node_num, false);
   1.149        LargeCostVector hyper_cost(_res_node_num);
   1.150 +      int relabel_cnt = 0;
   1.151        for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
   1.152                                          1 : _epsilon / _alpha )
   1.153        {
   1.154 @@ -1293,9 +1312,11 @@
   1.155               min_red_cost = hyper[n] ? -hyper_cost[n] :
   1.156                 std::numeric_limits<LargeCost>::max();
   1.157              for (int a = _first_out[n]; a != last_out; ++a) {
   1.158 -              rc = _cost[a] + pi_n - _pi[_target[a]];
   1.159 -              if (_res_cap[a] > 0 && rc < min_red_cost) {
   1.160 -                min_red_cost = rc;
   1.161 +              if (_res_cap[a] > 0) {
   1.162 +                rc = _cost[a] + pi_n - _pi[_target[a]];
   1.163 +                if (rc < min_red_cost) {
   1.164 +                  min_red_cost = rc;
   1.165 +                }
   1.166                }
   1.167              }
   1.168              _pi[n] -= min_red_cost + _epsilon;
   1.169 @@ -1313,11 +1334,11 @@
   1.170            }
   1.171  
   1.172            // Global update heuristic
   1.173 -          if (relabel_cnt >= next_update_limit) {
   1.174 +          if (relabel_cnt >= next_global_update_limit) {
   1.175              globalUpdate();
   1.176              for (int u = 0; u != _res_node_num; ++u)
   1.177                hyper[u] = false;
   1.178 -            next_update_limit += global_update_freq;
   1.179 +            next_global_update_limit += global_update_skip;
   1.180            }
   1.181          }
   1.182        }