COIN-OR::LEMON - Graph Library

Changeset 2635:23570e368afa in lemon-0.x for lemon/network_simplex.h


Ignore:
Timestamp:
06/01/09 17:35:27 (10 years ago)
Author:
Peter Kovacs
Branch:
default
Phase:
public
Convert:
svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3520
Message:

XTI data structure for NS - backport hg commit [8c3112a66878]

File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/network_simplex.h

    r2634 r2635  
    242242        _edge(ns._edge), _source(ns._source), _target(ns._target),
    243243        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
    244         _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
     244        _in_edge(ns._in_edge), _edge_num(ns._edge_num + ns._node_num), _next_edge(0)
    245245      {
    246246        // The main parameters of the pivot rule
     
    592592    int _edge_num;
    593593
    594     // The depth vector of the spanning tree structure
    595     IntVector _depth;
    596594    // The parent vector of the spanning tree structure
    597595    IntVector _parent;
     
    600598    // The thread vector of the spanning tree structure
    601599    IntVector _thread;
     600
     601    IntVector _rev_thread;
     602    IntVector _succ_num;
     603    IntVector _last_succ;
     604
     605    IntVector _dirty_revs;
     606   
    602607    // The forward vector of the spanning tree structure
    603608    BoolVector _forward;
     
    883888      _pi.resize(all_node_num, 0);
    884889     
    885       _depth.resize(all_node_num);
    886890      _parent.resize(all_node_num);
    887891      _pred.resize(all_node_num);
     892      _forward.resize(all_node_num);
    888893      _thread.resize(all_node_num);
    889       _forward.resize(all_node_num);
     894      _rev_thread.resize(all_node_num);
     895      _succ_num.resize(all_node_num);
     896      _last_succ.resize(all_node_num);
    890897      _state.resize(all_edge_num, STATE_LOWER);
    891898     
     
    916923      // Set data for the artificial root node
    917924      _root = _node_num;
    918       _depth[_root] = 0;
    919925      _parent[_root] = -1;
    920926      _pred[_root] = -1;
    921927      _thread[_root] = 0;
     928      _rev_thread[0] = _root;
     929      _succ_num[_root] = all_node_num;
     930      _last_succ[_root] = _root - 1;
    922931      _supply[_root] = 0;
    923932      _pi[_root] = 0;
     
    956965      Capacity max_cap = std::numeric_limits<Capacity>::max();
    957966      for (int u = 0, e = _edge_num; u != _node_num; ++u, ++e) {
    958         _thread[u] = u + 1;
    959         _depth[u] = 1;
    960967        _parent[u] = _root;
    961968        _pred[u] = e;
     969        _thread[u] = u + 1;
     970        _rev_thread[u + 1] = u;
     971        _succ_num[u] = 1;
     972        _last_succ[u] = u;
     973        _cap[e] = max_cap;
     974        _state[e] = STATE_TREE;
    962975        if (_supply[u] >= 0) {
     976          _forward[u] = true;
     977          _pi[u] = 0;
     978          _source[e] = u;
     979          _target[e] = _root;
    963980          _flow[e] = _supply[u];
    964           _forward[u] = true;
    965           _pi[u] = -max_cost;
    966         } else {
    967           _flow[e] = -_supply[u];
     981          _cost[e] = 0;
     982        }
     983        else {
    968984          _forward[u] = false;
    969985          _pi[u] = max_cost;
    970         }
    971         _cost[e] = max_cost;
    972         _cap[e] = max_cap;
    973         _state[e] = STATE_TREE;
     986          _source[e] = _root;
     987          _target[e] = u;
     988          _flow[e] = -_supply[u];
     989          _cost[e] = max_cost;
     990        }
    974991      }
    975992
     
    981998      int u = _source[_in_edge];
    982999      int v = _target[_in_edge];
    983       while (_depth[u] > _depth[v]) u = _parent[u];
    984       while (_depth[v] > _depth[u]) v = _parent[v];
    9851000      while (u != v) {
    986         u = _parent[u];
    987         v = _parent[v];
     1001        if (_succ_num[u] < _succ_num[v]) {
     1002          u = _parent[u];
     1003        } else {
     1004          v = _parent[v];
     1005        }
    9881006      }
    9891007      join = u;
     
    10601078      }
    10611079    }
    1062 
    1063     // Update _thread and _parent vectors
    1064     void updateThreadParent() {
    1065       int u;
     1080   
     1081    // Update the tree structure
     1082    void updateTreeStructure() {
     1083      int u, w;
     1084      int old_rev_thread = _rev_thread[u_out];
     1085      int old_succ_num = _succ_num[u_out];
     1086      int old_last_succ = _last_succ[u_out];
    10661087      v_out = _parent[u_out];
    10671088
    1068       // Handle the case when join and v_out coincide
    1069       bool par_first = false;
    1070       if (join == v_out) {
    1071         for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
    1072         if (u == v_in) {
    1073           par_first = true;
    1074           while (_thread[u] != u_out) u = _thread[u];
    1075           first = u;
    1076         }
    1077       }
    1078 
    1079       // Find the last successor of u_in (u) and the node after it (right)
    1080       // according to the thread index
    1081       for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
    1082       right = _thread[u];
    1083       if (_thread[v_in] == u_out) {
    1084         for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;
    1085         if (last == u_out) last = _thread[last];
    1086       }
    1087       else last = _thread[v_in];
    1088 
    1089       // Update stem nodes
     1089      u = _last_succ[u_in];  // the last successor of u_in
     1090      right = _thread[u];    // the node after it
     1091     
     1092      // Handle the case when old_rev_thread equals to v_in
     1093      // (it also means that join and v_out coincide)
     1094      if (old_rev_thread == v_in) {
     1095        last = _thread[_last_succ[u_out]];
     1096      } else {
     1097        last = _thread[v_in];
     1098      }
     1099
     1100      // Update _thread and _parent along the stem nodes (i.e. the nodes
     1101      // between u_in and u_out, whose parent have to be changed)
    10901102      _thread[v_in] = stem = u_in;
     1103      _dirty_revs.clear();
     1104      _dirty_revs.push_back(v_in);
    10911105      par_stem = v_in;
    10921106      while (stem != u_out) {
    1093         _thread[u] = new_stem = _parent[stem];
    1094 
    1095         // Find the node just before the stem node (u) according to
    1096         // the original thread index
    1097         for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
    1098         _thread[u] = right;
    1099 
    1100         // Change the parent node of stem and shift stem and par_stem nodes
    1101         _parent[stem] = par_stem;
     1107        // Insert the next stem node into the thread list
     1108        new_stem = _parent[stem];
     1109        _thread[u] = new_stem;
     1110        _dirty_revs.push_back(u);
     1111       
     1112        // Remove the subtree of stem from the thread list
     1113        w = _rev_thread[stem];
     1114        _thread[w] = right;
     1115        _rev_thread[right] = w;
     1116
     1117        // Change the parent node and shift stem nodes
     1118        _parent[stem] = par_stem;       
    11021119        par_stem = stem;
    11031120        stem = new_stem;
    11041121
    1105         // Find the last successor of stem (u) and the node after it (right)
    1106         // according to the thread index
    1107         for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
     1122        // Update u and right nodes
     1123        u = _last_succ[stem] == _last_succ[par_stem] ?
     1124          _rev_thread[par_stem] : _last_succ[stem];
    11081125        right = _thread[u];
    11091126      }
    11101127      _parent[u_out] = par_stem;
     1128      _last_succ[u_out] = u;
    11111129      _thread[u] = last;
    1112 
    1113       if (join == v_out && par_first) {
    1114         if (first != v_in) _thread[first] = right;
     1130      _rev_thread[last] = u;
     1131
     1132      // Remove the subtree of u_out from the thread list except for
     1133      // the case when old_rev_thread equals to v_in
     1134      // (it also means that join and v_out coincide)
     1135      if (old_rev_thread != v_in) {
     1136        _thread[old_rev_thread] = right;
     1137        _rev_thread[right] = old_rev_thread;
     1138      }
     1139     
     1140      // Update _rev_thread using the new _thread values
     1141      for (int i = 0; i < int(_dirty_revs.size()); ++i) {
     1142        u = _dirty_revs[i];
     1143        _rev_thread[_thread[u]] = u;
     1144      }
     1145     
     1146      // Update _last_succ for the stem nodes from u_out to u_in
     1147      for (u = u_out; u != u_in; u = _parent[u]) {
     1148        _last_succ[_parent[u]] = _last_succ[u];     
     1149      }
     1150
     1151      // Set limits for updating _last_succ form v_in and v_out
     1152      // towards the root
     1153      int up_limit_in = -1;
     1154      int up_limit_out = -1;
     1155      if (_last_succ[join] == v_in) {
     1156        up_limit_out = join;
    11151157      } else {
    1116         for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
    1117         _thread[u] = right;
    1118       }
    1119     }
    1120 
    1121     // Update _pred and _forward vectors
    1122     void updatePredEdge() {
    1123       int u = u_out, v;
     1158        up_limit_in = join;
     1159      }
     1160
     1161      // Update _last_succ from v_in towards the root
     1162      for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
     1163           u = _parent[u]) {
     1164        _last_succ[u] = _last_succ[u_out];
     1165      }
     1166      // Update _last_succ from v_out towards the root
     1167      if (join != old_rev_thread && v_in != old_rev_thread) {
     1168        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
     1169             u = _parent[u]) {
     1170          _last_succ[u] = old_rev_thread;
     1171        }
     1172      } else {
     1173        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
     1174             u = _parent[u]) {
     1175          _last_succ[u] = _last_succ[u_out];
     1176        }
     1177      }
     1178
     1179      // Update _pred and _forward for the stem nodes from u_out to u_in
     1180      u = u_out;
    11241181      while (u != u_in) {
    1125         v = _parent[u];
    1126         _pred[u] = _pred[v];
    1127         _forward[u] = !_forward[v];
    1128         u = v;
     1182        w = _parent[u];
     1183        _pred[u] = _pred[w];
     1184        _forward[u] = !_forward[w];
     1185        u = w;
    11291186      }
    11301187      _pred[u_in] = _in_edge;
    11311188      _forward[u_in] = (u_in == _source[_in_edge]);
    1132     }
    1133 
    1134     // Update _depth and _potential vectors
    1135     void updateDepthPotential() {
    1136       _depth[u_in] = _depth[v_in] + 1;
     1189     
     1190      // Update _succ_num from v_in to join
     1191      for (u = v_in; u != join; u = _parent[u]) {
     1192        _succ_num[u] += old_succ_num;
     1193      }
     1194      // Update _succ_num from v_out to join
     1195      for (u = v_out; u != join; u = _parent[u]) {
     1196        _succ_num[u] -= old_succ_num;
     1197      }
     1198      // Update _succ_num for the stem nodes from u_out to u_in
     1199      int tmp = 0;
     1200      u = u_out;
     1201      while (u != u_in) {
     1202        w = _parent[u];
     1203        tmp = _succ_num[u] - _succ_num[w] + tmp;
     1204        _succ_num[u] = tmp;
     1205        u = w;
     1206      }
     1207      _succ_num[u_in] = old_succ_num;
     1208    }
     1209
     1210    // Update potentials
     1211    void updatePotential() {
    11371212      Cost sigma = _forward[u_in] ?
    11381213        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
    11391214        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
    1140       _pi[u_in] += sigma;
    1141       for(int u = _thread[u_in]; _parent[u] != -1; u = _thread[u]) {
    1142         _depth[u] = _depth[_parent[u]] + 1;
    1143         if (_depth[u] <= _depth[u_in]) break;
     1215      // Update in the lower subtree (which has been moved)
     1216      int end = _thread[_last_succ[u_in]];
     1217      for (int u = u_in; u != end; u = _thread[u]) {
    11441218        _pi[u] += sigma;
    11451219      }
     
    11671241    bool start() {
    11681242      PivotRuleImplementation pivot(*this);
    1169 
     1243     
    11701244      // Execute the network simplex algorithm
    11711245      while (pivot.findEnteringEdge()) {
     
    11741248        changeFlow(change);
    11751249        if (change) {
    1176           updateThreadParent();
    1177           updatePredEdge();
    1178           updateDepthPotential();
    1179         }
    1180       }
    1181 
     1250          updateTreeStructure();
     1251          updatePotential();
     1252        }
     1253      }
     1254     
    11821255      // Check if the flow amount equals zero on all the artificial edges
    11831256      for (int e = _edge_num; e != _edge_num + _node_num; ++e) {
     
    12001273        (*_potential_result)[_node[i]] = _pi[i];
    12011274      }
    1202 
     1275     
    12031276      return true;
    12041277    }
Note: See TracChangeset for help on using the changeset viewer.