lemon/cost_scaling.h
changeset 936 ddd3c0d3d9bf
parent 935 6ea176638264
child 937 1226290a9b7d
     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