author Peter Kovacs Tue, 15 Mar 2011 19:32:21 +0100 changeset 1047 ddd3c0d3d9bf parent 1046 6ea176638264 child 1048 1226290a9b7d
Implement the scaling Price Refinement heuristic in CostScaling (#417)

These two heuristics are similar, but the newer one is faster
and not only makes it possible to skip some epsilon phases, but
it can improve the performance of the other phases, as well.
 lemon/cost_scaling.h file | annotate | diff | comparison | revisions
```     1.1 --- a/lemon/cost_scaling.h	Tue Mar 15 19:16:20 2011 +0100
1.2 +++ b/lemon/cost_scaling.h	Tue Mar 15 19:32:21 2011 +0100
1.3 @@ -980,30 +980,229 @@
1.4        }
1.5      }
1.6
1.7 -    // Early termination heuristic
1.8 -    bool earlyTermination() {
1.9 -      const double EARLY_TERM_FACTOR = 3.0;
1.10 +    // Price (potential) refinement heuristic
1.11 +    bool priceRefinement() {
1.12
1.13 -      // Build a static residual graph
1.14 -      _arc_vec.clear();
1.15 -      _cost_vec.clear();
1.16 -      for (int j = 0; j != _res_arc_num; ++j) {
1.17 -        if (_res_cap[j] > 0) {
1.18 -          _arc_vec.push_back(IntPair(_source[j], _target[j]));
1.19 -          _cost_vec.push_back(_cost[j] + 1);
1.20 +      // Stack for stroing the topological order
1.21 +      IntVector stack(_res_node_num);
1.22 +      int stack_top;
1.23 +
1.24 +      // Perform phases
1.25 +      while (topologicalSort(stack, stack_top)) {
1.26 +
1.27 +        // Compute node ranks in the acyclic admissible network and
1.28 +        // store the nodes in buckets
1.29 +        for (int i = 0; i != _res_node_num; ++i) {
1.30 +          _rank[i] = 0;
1.31          }
1.32 +        const int bucket_end = _root + 1;
1.33 +        for (int r = 0; r != _max_rank; ++r) {
1.34 +          _buckets[r] = bucket_end;
1.35 +        }
1.36 +        int top_rank = 0;
1.37 +        for ( ; stack_top >= 0; --stack_top) {
1.38 +          int u = stack[stack_top], v;
1.39 +          int rank_u = _rank[u];
1.40 +
1.41 +          LargeCost rc, pi_u = _pi[u];
1.42 +          int last_out = _first_out[u+1];
1.43 +          for (int a = _first_out[u]; a != last_out; ++a) {
1.44 +            if (_res_cap[a] > 0) {
1.45 +              v = _target[a];
1.46 +              rc = _cost[a] + pi_u - _pi[v];
1.47 +              if (rc < 0) {
1.48 +                LargeCost nrc = static_cast<LargeCost>((-rc - 0.5) / _epsilon);
1.49 +                if (nrc < LargeCost(_max_rank)) {
1.50 +                  int new_rank_v = rank_u + static_cast<int>(nrc);
1.51 +                  if (new_rank_v > _rank[v]) {
1.52 +                    _rank[v] = new_rank_v;
1.53 +                  }
1.54 +                }
1.55 +              }
1.56 +            }
1.57 +          }
1.58 +
1.59 +          if (rank_u > 0) {
1.60 +            top_rank = std::max(top_rank, rank_u);
1.61 +            int bfirst = _buckets[rank_u];
1.62 +            _bucket_next[u] = bfirst;
1.63 +            _bucket_prev[bfirst] = u;
1.64 +            _buckets[rank_u] = u;
1.65 +          }
1.66 +        }
1.67 +
1.68 +        // Check if the current flow is epsilon-optimal
1.69 +        if (top_rank == 0) {
1.70 +          return true;
1.71 +        }
1.72 +
1.73 +        // Process buckets in top-down order
1.74 +        for (int rank = top_rank; rank > 0; --rank) {
1.75 +          while (_buckets[rank] != bucket_end) {
1.76 +            // Remove the first node from the current bucket
1.77 +            int u = _buckets[rank];
1.78 +            _buckets[rank] = _bucket_next[u];
1.79 +
1.80 +            // Search the outgoing arcs of u
1.81 +            LargeCost rc, pi_u = _pi[u];
1.82 +            int last_out = _first_out[u+1];
1.83 +            int v, old_rank_v, new_rank_v;
1.84 +            for (int a = _first_out[u]; a != last_out; ++a) {
1.85 +              if (_res_cap[a] > 0) {
1.86 +                v = _target[a];
1.87 +                old_rank_v = _rank[v];
1.88 +
1.89 +                if (old_rank_v < rank) {
1.90 +
1.91 +                  // Compute the new rank of node v
1.92 +                  rc = _cost[a] + pi_u - _pi[v];
1.93 +                  if (rc < 0) {
1.94 +                    new_rank_v = rank;
1.95 +                  } else {
1.96 +                    LargeCost nrc = rc / _epsilon;
1.97 +                    new_rank_v = 0;
1.98 +                    if (nrc < LargeCost(_max_rank)) {
1.99 +                      new_rank_v = rank - 1 - static_cast<int>(nrc);
1.100 +                    }
1.101 +                  }
1.102 +
1.103 +                  // Change the rank of node v
1.104 +                  if (new_rank_v > old_rank_v) {
1.105 +                    _rank[v] = new_rank_v;
1.106 +
1.107 +                    // Remove v from its old bucket
1.108 +                    if (old_rank_v > 0) {
1.109 +                      if (_buckets[old_rank_v] == v) {
1.110 +                        _buckets[old_rank_v] = _bucket_next[v];
1.111 +                      } else {
1.112 +                        int pv = _bucket_prev[v], nv = _bucket_next[v];
1.113 +                        _bucket_next[pv] = nv;
1.114 +                        _bucket_prev[nv] = pv;
1.115 +                      }
1.116 +                    }
1.117 +
1.118 +                    // Insert v into its new bucket
1.119 +                    int nv = _buckets[new_rank_v];
1.120 +                    _bucket_next[v] = nv;
1.121 +                    _bucket_prev[nv] = v;
1.122 +                    _buckets[new_rank_v] = v;
1.123 +                  }
1.124 +                }
1.125 +              }
1.126 +            }
1.127 +
1.128 +            // Refine potential of node u
1.129 +            _pi[u] -= rank * _epsilon;
1.130 +          }
1.131 +        }
1.132 +
1.133        }
1.134 -      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1.135
1.136 -      // Run Bellman-Ford algorithm to check if the current flow is optimal
1.137 -      BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
1.138 -      bf.init(0);
1.139 -      bool done = false;
1.140 -      int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num)));
1.141 -      for (int i = 0; i < K && !done; ++i) {
1.142 -        done = bf.processNextWeakRound();
1.143 +      return false;
1.144 +    }
1.145 +
1.146 +    // Find and cancel cycles in the admissible network and
1.147 +    // determine topological order using DFS
1.148 +    bool topologicalSort(IntVector &stack, int &stack_top) {
1.149 +      const int MAX_CYCLE_CANCEL = 1;
1.150 +
1.151 +      BoolVector reached(_res_node_num, false);
1.152 +      BoolVector processed(_res_node_num, false);
1.153 +      IntVector pred(_res_node_num);
1.154 +      for (int i = 0; i != _res_node_num; ++i) {
1.155 +        _next_out[i] = _first_out[i];
1.156        }
1.157 -      return done;
1.158 +      stack_top = -1;
1.159 +
1.160 +      int cycle_cnt = 0;
1.161 +      for (int start = 0; start != _res_node_num; ++start) {
1.162 +        if (reached[start]) continue;
1.163 +
1.164 +        // Start DFS search from this start node
1.165 +        pred[start] = -1;
1.166 +        int tip = start, v;
1.167 +        while (true) {
1.168 +          // Check the outgoing arcs of the current tip node
1.169 +          reached[tip] = true;
1.170 +          LargeCost pi_tip = _pi[tip];
1.171 +          int a, last_out = _first_out[tip+1];
1.172 +          for (a = _next_out[tip]; a != last_out; ++a) {
1.173 +            if (_res_cap[a] > 0) {
1.174 +              v = _target[a];
1.175 +              if (_cost[a] + pi_tip - _pi[v] < 0) {
1.176 +                if (!reached[v]) {
1.177 +                  // A new node is reached
1.178 +                  reached[v] = true;
1.179 +                  pred[v] = tip;
1.180 +                  _next_out[tip] = a;
1.181 +                  tip = v;
1.182 +                  a = _next_out[tip];
1.183 +                  last_out = _first_out[tip+1];
1.184 +                  break;
1.185 +                }
1.186 +                else if (!processed[v]) {
1.187 +                  // A cycle is found
1.188 +                  ++cycle_cnt;
1.189 +                  _next_out[tip] = a;
1.190 +
1.191 +                  // Find the minimum residual capacity along the cycle
1.192 +                  Value d, delta = _res_cap[a];
1.193 +                  int u, delta_node = tip;
1.194 +                  for (u = tip; u != v; ) {
1.195 +                    u = pred[u];
1.196 +                    d = _res_cap[_next_out[u]];
1.197 +                    if (d <= delta) {
1.198 +                      delta = d;
1.199 +                      delta_node = u;
1.200 +                    }
1.201 +                  }
1.202 +
1.203 +                  // Augment along the cycle
1.204 +                  _res_cap[a] -= delta;
1.205 +                  _res_cap[_reverse[a]] += delta;
1.206 +                  for (u = tip; u != v; ) {
1.207 +                    u = pred[u];
1.208 +                    int ca = _next_out[u];
1.209 +                    _res_cap[ca] -= delta;
1.210 +                    _res_cap[_reverse[ca]] += delta;
1.211 +                  }
1.212 +
1.213 +                  // Check the maximum number of cycle canceling
1.214 +                  if (cycle_cnt >= MAX_CYCLE_CANCEL) {
1.215 +                    return false;
1.216 +                  }
1.217 +
1.218 +                  // Roll back search to delta_node
1.219 +                  if (delta_node != tip) {
1.220 +                    for (u = tip; u != delta_node; u = pred[u]) {
1.221 +                      reached[u] = false;
1.222 +                    }
1.223 +                    tip = delta_node;
1.224 +                    a = _next_out[tip] + 1;
1.225 +                    last_out = _first_out[tip+1];
1.226 +                    break;
1.227 +                  }
1.228 +                }
1.229 +              }
1.230 +            }
1.231 +          }
1.232 +
1.233 +          // Step back to the previous node
1.234 +          if (a == last_out) {
1.235 +            processed[tip] = true;
1.236 +            stack[++stack_top] = tip;
1.237 +            tip = pred[tip];
1.238 +            if (tip < 0) {
1.239 +              // Finish DFS from the current start node
1.240 +              break;
1.241 +            }
1.242 +            ++_next_out[tip];
1.243 +          }
1.244 +        }
1.245 +
1.246 +      }
1.247 +
1.248 +      return (cycle_cnt == 0);
1.249      }
1.250
1.251      // Global potential update heuristic
1.252 @@ -1102,7 +1301,7 @@
1.253      /// Execute the algorithm performing augment and relabel operations
1.254      void startAugment(int max_length) {
1.255        // Paramters for heuristics
1.256 -      const int EARLY_TERM_EPSILON_LIMIT = 1000;
1.257 +      const int PRICE_REFINEMENT_LIMIT = 2;
1.258        const double GLOBAL_UPDATE_FACTOR = 1.0;
1.259        const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
1.260          (_res_node_num + _sup_node_num * _sup_node_num));
1.261 @@ -1112,12 +1311,15 @@
1.262        IntVector path;
1.263        BoolVector path_arc(_res_arc_num, false);
1.264        int relabel_cnt = 0;
1.265 +      int eps_phase_cnt = 0;
1.266        for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.267                                          1 : _epsilon / _alpha )
1.268        {
1.269 -        // Early termination heuristic
1.270 -        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
1.271 -          if (earlyTermination()) break;
1.272 +        ++eps_phase_cnt;
1.273 +
1.274 +        // Price refinement heuristic
1.275 +        if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) {
1.276 +          if (priceRefinement()) continue;
1.277          }
1.278
1.279          // Initialize current phase
1.280 @@ -1224,7 +1426,7 @@
1.281      /// Execute the algorithm performing push and relabel operations
1.282      void startPush() {
1.283        // Paramters for heuristics
1.284 -      const int EARLY_TERM_EPSILON_LIMIT = 1000;
1.285 +      const int PRICE_REFINEMENT_LIMIT = 2;
1.286        const double GLOBAL_UPDATE_FACTOR = 2.0;
1.287
1.288        const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
1.289 @@ -1235,12 +1437,15 @@
1.290        BoolVector hyper(_res_node_num, false);
1.291        LargeCostVector hyper_cost(_res_node_num);
1.292        int relabel_cnt = 0;
1.293 +      int eps_phase_cnt = 0;
1.294        for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.295                                          1 : _epsilon / _alpha )
1.296        {
1.297 -        // Early termination heuristic
1.298 -        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
1.299 -          if (earlyTermination()) break;
1.300 +        ++eps_phase_cnt;
1.301 +
1.302 +        // Price refinement heuristic
1.303 +        if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) {
1.304 +          if (priceRefinement()) continue;
1.305          }
1.306
1.307          // Initialize current phase
```