Improve the tree update process and a pivot rule (#391)
authorPeter Kovacs <kpeter@inf.elte.hu>
Sun, 22 Aug 2010 23:54:10 +0200
changeset 990dca9eed2c375
parent 989 24b3f18ed9e2
child 991 fb932bcfd803
Improve the tree update process and a pivot rule (#391)
and make some parts of the code clearer using better names
lemon/network_simplex.h
     1.1 --- a/lemon/network_simplex.h	Tue Jun 22 16:13:00 2010 +0200
     1.2 +++ b/lemon/network_simplex.h	Sun Aug 22 23:54:10 2010 +0200
     1.3 @@ -166,8 +166,9 @@
     1.4      typedef std::vector<int> IntVector;
     1.5      typedef std::vector<Value> ValueVector;
     1.6      typedef std::vector<Cost> CostVector;
     1.7 -    typedef std::vector<char> BoolVector;
     1.8 -    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
     1.9 +    typedef std::vector<signed char> CharVector;
    1.10 +    // Note: vector<signed char> is used instead of vector<ArcState> and 
    1.11 +    // vector<ArcDirection> for efficiency reasons
    1.12  
    1.13      // State constants for arcs
    1.14      enum ArcState {
    1.15 @@ -176,9 +177,11 @@
    1.16        STATE_LOWER =  1
    1.17      };
    1.18  
    1.19 -    typedef std::vector<signed char> StateVector;
    1.20 -    // Note: vector<signed char> is used instead of vector<ArcState> for
    1.21 -    // efficiency reasons
    1.22 +    // Direction constants for tree arcs
    1.23 +    enum ArcDirection {
    1.24 +      DIR_DOWN = -1,
    1.25 +      DIR_UP   =  1
    1.26 +    };
    1.27  
    1.28    private:
    1.29  
    1.30 @@ -217,15 +220,13 @@
    1.31      IntVector _rev_thread;
    1.32      IntVector _succ_num;
    1.33      IntVector _last_succ;
    1.34 +    CharVector _pred_dir;
    1.35 +    CharVector _state;
    1.36      IntVector _dirty_revs;
    1.37 -    BoolVector _forward;
    1.38 -    StateVector _state;
    1.39      int _root;
    1.40  
    1.41      // Temporary data used in the current pivot iteration
    1.42      int in_arc, join, u_in, v_in, u_out, v_out;
    1.43 -    int first, second, right, last;
    1.44 -    int stem, par_stem, new_stem;
    1.45      Value delta;
    1.46  
    1.47      const Value MAX;
    1.48 @@ -250,7 +251,7 @@
    1.49        const IntVector  &_source;
    1.50        const IntVector  &_target;
    1.51        const CostVector &_cost;
    1.52 -      const StateVector &_state;
    1.53 +      const CharVector &_state;
    1.54        const CostVector &_pi;
    1.55        int &_in_arc;
    1.56        int _search_arc_num;
    1.57 @@ -302,7 +303,7 @@
    1.58        const IntVector  &_source;
    1.59        const IntVector  &_target;
    1.60        const CostVector &_cost;
    1.61 -      const StateVector &_state;
    1.62 +      const CharVector &_state;
    1.63        const CostVector &_pi;
    1.64        int &_in_arc;
    1.65        int _search_arc_num;
    1.66 @@ -341,7 +342,7 @@
    1.67        const IntVector  &_source;
    1.68        const IntVector  &_target;
    1.69        const CostVector &_cost;
    1.70 -      const StateVector &_state;
    1.71 +      const CharVector &_state;
    1.72        const CostVector &_pi;
    1.73        int &_in_arc;
    1.74        int _search_arc_num;
    1.75 @@ -414,7 +415,7 @@
    1.76        const IntVector  &_source;
    1.77        const IntVector  &_target;
    1.78        const CostVector &_cost;
    1.79 -      const StateVector &_state;
    1.80 +      const CharVector &_state;
    1.81        const CostVector &_pi;
    1.82        int &_in_arc;
    1.83        int _search_arc_num;
    1.84 @@ -517,7 +518,7 @@
    1.85        const IntVector  &_source;
    1.86        const IntVector  &_target;
    1.87        const CostVector &_cost;
    1.88 -      const StateVector &_state;
    1.89 +      const CharVector &_state;
    1.90        const CostVector &_pi;
    1.91        int &_in_arc;
    1.92        int _search_arc_num;
    1.93 @@ -570,11 +571,13 @@
    1.94        bool findEnteringArc() {
    1.95          // Check the current candidate list
    1.96          int e;
    1.97 +        Cost c;
    1.98          for (int i = 0; i != _curr_length; ++i) {
    1.99            e = _candidates[i];
   1.100 -          _cand_cost[e] = _state[e] *
   1.101 -            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   1.102 -          if (_cand_cost[e] >= 0) {
   1.103 +          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   1.104 +          if (c < 0) {
   1.105 +            _cand_cost[e] = c;
   1.106 +          } else {
   1.107              _candidates[i--] = _candidates[--_curr_length];
   1.108            }
   1.109          }
   1.110 @@ -584,9 +587,9 @@
   1.111          int limit = _head_length;
   1.112  
   1.113          for (e = _next_arc; e != _search_arc_num; ++e) {
   1.114 -          _cand_cost[e] = _state[e] *
   1.115 -            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   1.116 -          if (_cand_cost[e] < 0) {
   1.117 +          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   1.118 +          if (c < 0) {
   1.119 +            _cand_cost[e] = c;
   1.120              _candidates[_curr_length++] = e;
   1.121            }
   1.122            if (--cnt == 0) {
   1.123 @@ -913,7 +916,7 @@
   1.124  
   1.125        _parent.resize(all_node_num);
   1.126        _pred.resize(all_node_num);
   1.127 -      _forward.resize(all_node_num);
   1.128 +      _pred_dir.resize(all_node_num);
   1.129        _thread.resize(all_node_num);
   1.130        _rev_thread.resize(all_node_num);
   1.131        _succ_num.resize(all_node_num);
   1.132 @@ -1116,14 +1119,14 @@
   1.133            _cap[e] = INF;
   1.134            _state[e] = STATE_TREE;
   1.135            if (_supply[u] >= 0) {
   1.136 -            _forward[u] = true;
   1.137 +            _pred_dir[u] = DIR_UP;
   1.138              _pi[u] = 0;
   1.139              _source[e] = u;
   1.140              _target[e] = _root;
   1.141              _flow[e] = _supply[u];
   1.142              _cost[e] = 0;
   1.143            } else {
   1.144 -            _forward[u] = false;
   1.145 +            _pred_dir[u] = DIR_DOWN;
   1.146              _pi[u] = ART_COST;
   1.147              _source[e] = _root;
   1.148              _target[e] = u;
   1.149 @@ -1143,7 +1146,7 @@
   1.150            _succ_num[u] = 1;
   1.151            _last_succ[u] = u;
   1.152            if (_supply[u] >= 0) {
   1.153 -            _forward[u] = true;
   1.154 +            _pred_dir[u] = DIR_UP;
   1.155              _pi[u] = 0;
   1.156              _pred[u] = e;
   1.157              _source[e] = u;
   1.158 @@ -1153,7 +1156,7 @@
   1.159              _cost[e] = 0;
   1.160              _state[e] = STATE_TREE;
   1.161            } else {
   1.162 -            _forward[u] = false;
   1.163 +            _pred_dir[u] = DIR_DOWN;
   1.164              _pi[u] = ART_COST;
   1.165              _pred[u] = f;
   1.166              _source[f] = _root;
   1.167 @@ -1184,7 +1187,7 @@
   1.168            _succ_num[u] = 1;
   1.169            _last_succ[u] = u;
   1.170            if (_supply[u] <= 0) {
   1.171 -            _forward[u] = false;
   1.172 +            _pred_dir[u] = DIR_DOWN;
   1.173              _pi[u] = 0;
   1.174              _pred[u] = e;
   1.175              _source[e] = _root;
   1.176 @@ -1194,7 +1197,7 @@
   1.177              _cost[e] = 0;
   1.178              _state[e] = STATE_TREE;
   1.179            } else {
   1.180 -            _forward[u] = true;
   1.181 +            _pred_dir[u] = DIR_UP;
   1.182              _pi[u] = -ART_COST;
   1.183              _pred[u] = f;
   1.184              _source[f] = u;
   1.185 @@ -1237,6 +1240,7 @@
   1.186      bool findLeavingArc() {
   1.187        // Initialize first and second nodes according to the direction
   1.188        // of the cycle
   1.189 +      int first, second;
   1.190        if (_state[in_arc] == STATE_LOWER) {
   1.191          first  = _source[in_arc];
   1.192          second = _target[in_arc];
   1.193 @@ -1246,25 +1250,32 @@
   1.194        }
   1.195        delta = _cap[in_arc];
   1.196        int result = 0;
   1.197 -      Value d;
   1.198 +      Value c, d;
   1.199        int e;
   1.200  
   1.201 -      // Search the cycle along the path form the first node to the root
   1.202 +      // Search the cycle form the first node to the join node
   1.203        for (int u = first; u != join; u = _parent[u]) {
   1.204          e = _pred[u];
   1.205 -        d = _forward[u] ?
   1.206 -          _flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]);
   1.207 +        d = _flow[e];
   1.208 +        if (_pred_dir[u] == DIR_DOWN) {
   1.209 +          c = _cap[e];
   1.210 +          d = c >= MAX ? INF : c - d;
   1.211 +        }
   1.212          if (d < delta) {
   1.213            delta = d;
   1.214            u_out = u;
   1.215            result = 1;
   1.216          }
   1.217        }
   1.218 -      // Search the cycle along the path form the second node to the root
   1.219 +
   1.220 +      // Search the cycle form the second node to the join node
   1.221        for (int u = second; u != join; u = _parent[u]) {
   1.222          e = _pred[u];
   1.223 -        d = _forward[u] ?
   1.224 -          (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
   1.225 +        d = _flow[e];
   1.226 +        if (_pred_dir[u] == DIR_UP) {
   1.227 +          c = _cap[e];
   1.228 +          d = c >= MAX ? INF : c - d;
   1.229 +        }
   1.230          if (d <= delta) {
   1.231            delta = d;
   1.232            u_out = u;
   1.233 @@ -1289,10 +1300,10 @@
   1.234          Value val = _state[in_arc] * delta;
   1.235          _flow[in_arc] += val;
   1.236          for (int u = _source[in_arc]; u != join; u = _parent[u]) {
   1.237 -          _flow[_pred[u]] += _forward[u] ? -val : val;
   1.238 +          _flow[_pred[u]] -= _pred_dir[u] * val;
   1.239          }
   1.240          for (int u = _target[in_arc]; u != join; u = _parent[u]) {
   1.241 -          _flow[_pred[u]] += _forward[u] ? val : -val;
   1.242 +          _flow[_pred[u]] += _pred_dir[u] * val;
   1.243          }
   1.244        }
   1.245        // Update the state of the entering and leaving arcs
   1.246 @@ -1307,130 +1318,134 @@
   1.247  
   1.248      // Update the tree structure
   1.249      void updateTreeStructure() {
   1.250 -      int u, w;
   1.251        int old_rev_thread = _rev_thread[u_out];
   1.252        int old_succ_num = _succ_num[u_out];
   1.253        int old_last_succ = _last_succ[u_out];
   1.254        v_out = _parent[u_out];
   1.255  
   1.256 -      u = _last_succ[u_in];  // the last successor of u_in
   1.257 -      right = _thread[u];    // the node after it
   1.258 +      // Check if u_in and u_out coincide
   1.259 +      if (u_in == u_out) {
   1.260 +        // Update _parent, _pred, _pred_dir
   1.261 +        _parent[u_in] = v_in;
   1.262 +        _pred[u_in] = in_arc;
   1.263 +        _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
   1.264  
   1.265 -      // Handle the case when old_rev_thread equals to v_in
   1.266 -      // (it also means that join and v_out coincide)
   1.267 -      if (old_rev_thread == v_in) {
   1.268 -        last = _thread[_last_succ[u_out]];
   1.269 +        // Update _thread and _rev_thread
   1.270 +        if (_thread[v_in] != u_out) {
   1.271 +          int after = _thread[old_last_succ];
   1.272 +          _thread[old_rev_thread] = after;
   1.273 +          _rev_thread[after] = old_rev_thread;
   1.274 +          after = _thread[v_in];
   1.275 +          _thread[v_in] = u_out;
   1.276 +          _rev_thread[u_out] = v_in;
   1.277 +          _thread[old_last_succ] = after;
   1.278 +          _rev_thread[after] = old_last_succ;
   1.279 +        }
   1.280        } else {
   1.281 -        last = _thread[v_in];
   1.282 -      }
   1.283 +        // Handle the case when old_rev_thread equals to v_in
   1.284 +        // (it also means that join and v_out coincide)
   1.285 +        int thread_continue = old_rev_thread == v_in ?
   1.286 +          _thread[old_last_succ] : _thread[v_in];
   1.287  
   1.288 -      // Update _thread and _parent along the stem nodes (i.e. the nodes
   1.289 -      // between u_in and u_out, whose parent have to be changed)
   1.290 -      _thread[v_in] = stem = u_in;
   1.291 -      _dirty_revs.clear();
   1.292 -      _dirty_revs.push_back(v_in);
   1.293 -      par_stem = v_in;
   1.294 -      while (stem != u_out) {
   1.295 -        // Insert the next stem node into the thread list
   1.296 -        new_stem = _parent[stem];
   1.297 -        _thread[u] = new_stem;
   1.298 -        _dirty_revs.push_back(u);
   1.299 +        // Update _thread and _parent along the stem nodes (i.e. the nodes
   1.300 +        // between u_in and u_out, whose parent have to be changed)
   1.301 +        int stem = u_in;              // the current stem node
   1.302 +        int par_stem = v_in;          // the new parent of stem
   1.303 +        int next_stem;                // the next stem node
   1.304 +        int last = _last_succ[u_in];  // the last successor of stem
   1.305 +        int before, after = _thread[last];
   1.306 +        _thread[v_in] = u_in;
   1.307 +        _dirty_revs.clear();
   1.308 +        _dirty_revs.push_back(v_in);
   1.309 +        while (stem != u_out) {
   1.310 +          // Insert the next stem node into the thread list
   1.311 +          next_stem = _parent[stem];
   1.312 +          _thread[last] = next_stem;
   1.313 +          _dirty_revs.push_back(last);
   1.314  
   1.315 -        // Remove the subtree of stem from the thread list
   1.316 -        w = _rev_thread[stem];
   1.317 -        _thread[w] = right;
   1.318 -        _rev_thread[right] = w;
   1.319 +          // Remove the subtree of stem from the thread list
   1.320 +          before = _rev_thread[stem];
   1.321 +          _thread[before] = after;
   1.322 +          _rev_thread[after] = before;
   1.323  
   1.324 -        // Change the parent node and shift stem nodes
   1.325 -        _parent[stem] = par_stem;
   1.326 -        par_stem = stem;
   1.327 -        stem = new_stem;
   1.328 +          // Change the parent node and shift stem nodes
   1.329 +          _parent[stem] = par_stem;
   1.330 +          par_stem = stem;
   1.331 +          stem = next_stem;
   1.332  
   1.333 -        // Update u and right
   1.334 -        u = _last_succ[stem] == _last_succ[par_stem] ?
   1.335 -          _rev_thread[par_stem] : _last_succ[stem];
   1.336 -        right = _thread[u];
   1.337 -      }
   1.338 -      _parent[u_out] = par_stem;
   1.339 -      _thread[u] = last;
   1.340 -      _rev_thread[last] = u;
   1.341 -      _last_succ[u_out] = u;
   1.342 +          // Update last and after
   1.343 +          last = _last_succ[stem] == _last_succ[par_stem] ?
   1.344 +            _rev_thread[par_stem] : _last_succ[stem];
   1.345 +          after = _thread[last];
   1.346 +        }
   1.347 +        _parent[u_out] = par_stem;
   1.348 +        _thread[last] = thread_continue;
   1.349 +        _rev_thread[thread_continue] = last;
   1.350 +        _last_succ[u_out] = last;
   1.351  
   1.352 -      // Remove the subtree of u_out from the thread list except for
   1.353 -      // the case when old_rev_thread equals to v_in
   1.354 -      // (it also means that join and v_out coincide)
   1.355 -      if (old_rev_thread != v_in) {
   1.356 -        _thread[old_rev_thread] = right;
   1.357 -        _rev_thread[right] = old_rev_thread;
   1.358 -      }
   1.359 +        // Remove the subtree of u_out from the thread list except for
   1.360 +        // the case when old_rev_thread equals to v_in
   1.361 +        if (old_rev_thread != v_in) {
   1.362 +          _thread[old_rev_thread] = after;
   1.363 +          _rev_thread[after] = old_rev_thread;
   1.364 +        }
   1.365  
   1.366 -      // Update _rev_thread using the new _thread values
   1.367 -      for (int i = 0; i != int(_dirty_revs.size()); ++i) {
   1.368 -        u = _dirty_revs[i];
   1.369 -        _rev_thread[_thread[u]] = u;
   1.370 -      }
   1.371 +        // Update _rev_thread using the new _thread values
   1.372 +        for (int i = 0; i != int(_dirty_revs.size()); ++i) {
   1.373 +          int u = _dirty_revs[i];
   1.374 +          _rev_thread[_thread[u]] = u;
   1.375 +        }
   1.376  
   1.377 -      // Update _pred, _forward, _last_succ and _succ_num for the
   1.378 -      // stem nodes from u_out to u_in
   1.379 -      int tmp_sc = 0, tmp_ls = _last_succ[u_out];
   1.380 -      u = u_out;
   1.381 -      while (u != u_in) {
   1.382 -        w = _parent[u];
   1.383 -        _pred[u] = _pred[w];
   1.384 -        _forward[u] = !_forward[w];
   1.385 -        tmp_sc += _succ_num[u] - _succ_num[w];
   1.386 -        _succ_num[u] = tmp_sc;
   1.387 -        _last_succ[w] = tmp_ls;
   1.388 -        u = w;
   1.389 -      }
   1.390 -      _pred[u_in] = in_arc;
   1.391 -      _forward[u_in] = (u_in == _source[in_arc]);
   1.392 -      _succ_num[u_in] = old_succ_num;
   1.393 -
   1.394 -      // Set limits for updating _last_succ form v_in and v_out
   1.395 -      // towards the root
   1.396 -      int up_limit_in = -1;
   1.397 -      int up_limit_out = -1;
   1.398 -      if (_last_succ[join] == v_in) {
   1.399 -        up_limit_out = join;
   1.400 -      } else {
   1.401 -        up_limit_in = join;
   1.402 +        // Update _pred, _pred_dir, _last_succ and _succ_num for the
   1.403 +        // stem nodes from u_out to u_in
   1.404 +        int tmp_sc = 0, tmp_ls = _last_succ[u_out];
   1.405 +        for (int u = u_out, p = _parent[u]; u != u_in; u = p, p = _parent[u]) {
   1.406 +          _pred[u] = _pred[p];
   1.407 +          _pred_dir[u] = -_pred_dir[p];
   1.408 +          tmp_sc += _succ_num[u] - _succ_num[p];
   1.409 +          _succ_num[u] = tmp_sc;
   1.410 +          _last_succ[p] = tmp_ls;
   1.411 +        }
   1.412 +        _pred[u_in] = in_arc;
   1.413 +        _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
   1.414 +        _succ_num[u_in] = old_succ_num;
   1.415        }
   1.416  
   1.417        // Update _last_succ from v_in towards the root
   1.418 -      for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
   1.419 -           u = _parent[u]) {
   1.420 -        _last_succ[u] = _last_succ[u_out];
   1.421 +      int up_limit_out = _last_succ[join] == v_in ? join : -1;
   1.422 +      int last_succ_out = _last_succ[u_out];
   1.423 +      for (int u = v_in; u != -1 && _last_succ[u] == v_in; u = _parent[u]) {
   1.424 +        _last_succ[u] = last_succ_out;
   1.425        }
   1.426 +
   1.427        // Update _last_succ from v_out towards the root
   1.428        if (join != old_rev_thread && v_in != old_rev_thread) {
   1.429 -        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
   1.430 +        for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
   1.431               u = _parent[u]) {
   1.432            _last_succ[u] = old_rev_thread;
   1.433          }
   1.434 -      } else {
   1.435 -        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
   1.436 +      }
   1.437 +      else if (last_succ_out != old_last_succ) {
   1.438 +        for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
   1.439               u = _parent[u]) {
   1.440 -          _last_succ[u] = _last_succ[u_out];
   1.441 +          _last_succ[u] = last_succ_out;
   1.442          }
   1.443        }
   1.444  
   1.445        // Update _succ_num from v_in to join
   1.446 -      for (u = v_in; u != join; u = _parent[u]) {
   1.447 +      for (int u = v_in; u != join; u = _parent[u]) {
   1.448          _succ_num[u] += old_succ_num;
   1.449        }
   1.450        // Update _succ_num from v_out to join
   1.451 -      for (u = v_out; u != join; u = _parent[u]) {
   1.452 +      for (int u = v_out; u != join; u = _parent[u]) {
   1.453          _succ_num[u] -= old_succ_num;
   1.454        }
   1.455      }
   1.456  
   1.457 -    // Update potentials
   1.458 +    // Update potentials in the subtree that has been moved
   1.459      void updatePotential() {
   1.460 -      Cost sigma = _forward[u_in] ?
   1.461 -        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
   1.462 -        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
   1.463 -      // Update potentials in the subtree, which has been moved
   1.464 +      Cost sigma = _pi[v_in] - _pi[u_in] -
   1.465 +                   _pred_dir[u_in] * _cost[in_arc];
   1.466        int end = _thread[_last_succ[u_in]];
   1.467        for (int u = u_in; u != end; u = _thread[u]) {
   1.468          _pi[u] += sigma;