COIN-OR::LEMON - Graph Library

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/network_simplex.h

    r978 r1026  
    4848  /// flow problem.
    4949  ///
    50   /// In general, %NetworkSimplex is the fastest implementation available
    51   /// in LEMON for this problem.
    52   /// Moreover, it supports both directions of the supply/demand inequality
    53   /// constraints. For more information, see \ref SupplyType.
     50  /// In general, \ref NetworkSimplex and \ref CostScaling are the fastest
     51  /// implementations available in LEMON for this problem.
     52  /// Furthermore, this class supports both directions of the supply/demand
     53  /// inequality constraints. For more information, see \ref SupplyType.
    5454  ///
    5555  /// Most of the parameters of the problem (except for the digraph)
     
    6464  /// algorithm. By default, it is the same as \c V.
    6565  ///
    66   /// \warning Both number types must be signed and all input data must
     66  /// \warning Both \c V and \c C must be signed number types.
     67  /// \warning All input data (capacities, supply values, and costs) must
    6768  /// be integer.
    6869  ///
     
    126127    /// of the algorithm.
    127128    /// By default, \ref BLOCK_SEARCH "Block Search" is used, which
    128     /// proved to be the most efficient and the most robust on various
     129    /// turend out to be the most efficient and the most robust on various
    129130    /// test inputs.
    130131    /// However, another pivot rule can be selected using the \ref run()
     
    167168    typedef std::vector<Value> ValueVector;
    168169    typedef std::vector<Cost> CostVector;
    169     typedef std::vector<char> BoolVector;
    170     // Note: vector<char> is used instead of vector<bool> for efficiency reasons
     170    typedef std::vector<signed char> CharVector;
     171    // Note: vector<signed char> is used instead of vector<ArcState> and
     172    // vector<ArcDirection> for efficiency reasons
    171173
    172174    // State constants for arcs
     
    177179    };
    178180
    179     typedef std::vector<signed char> StateVector;
    180     // Note: vector<signed char> is used instead of vector<ArcState> for
    181     // efficiency reasons
     181    // Direction constants for tree arcs
     182    enum ArcDirection {
     183      DIR_DOWN = -1,
     184      DIR_UP   =  1
     185    };
    182186
    183187  private:
     
    218222    IntVector _succ_num;
    219223    IntVector _last_succ;
     224    CharVector _pred_dir;
     225    CharVector _state;
    220226    IntVector _dirty_revs;
    221     BoolVector _forward;
    222     StateVector _state;
    223227    int _root;
    224228
    225229    // Temporary data used in the current pivot iteration
    226230    int in_arc, join, u_in, v_in, u_out, v_out;
    227     int first, second, right, last;
    228     int stem, par_stem, new_stem;
    229231    Value delta;
    230232
     
    251253      const IntVector  &_target;
    252254      const CostVector &_cost;
    253       const StateVector &_state;
     255      const CharVector &_state;
    254256      const CostVector &_pi;
    255257      int &_in_arc;
     
    303305      const IntVector  &_target;
    304306      const CostVector &_cost;
    305       const StateVector &_state;
     307      const CharVector &_state;
    306308      const CostVector &_pi;
    307309      int &_in_arc;
     
    342344      const IntVector  &_target;
    343345      const CostVector &_cost;
    344       const StateVector &_state;
     346      const CharVector &_state;
    345347      const CostVector &_pi;
    346348      int &_in_arc;
     
    415417      const IntVector  &_target;
    416418      const CostVector &_cost;
    417       const StateVector &_state;
     419      const CharVector &_state;
    418420      const CostVector &_pi;
    419421      int &_in_arc;
     
    518520      const IntVector  &_target;
    519521      const CostVector &_cost;
    520       const StateVector &_state;
     522      const CharVector &_state;
    521523      const CostVector &_pi;
    522524      int &_in_arc;
     
    571573        // Check the current candidate list
    572574        int e;
     575        Cost c;
    573576        for (int i = 0; i != _curr_length; ++i) {
    574577          e = _candidates[i];
    575           _cand_cost[e] = _state[e] *
    576             (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    577           if (_cand_cost[e] >= 0) {
     578          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     579          if (c < 0) {
     580            _cand_cost[e] = c;
     581          } else {
    578582            _candidates[i--] = _candidates[--_curr_length];
    579583          }
     
    585589
    586590        for (e = _next_arc; e != _search_arc_num; ++e) {
    587           _cand_cost[e] = _state[e] *
    588             (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    589           if (_cand_cost[e] < 0) {
     591          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     592          if (c < 0) {
     593            _cand_cost[e] = c;
    590594            _candidates[_curr_length++] = e;
    591595          }
     
    634638    ///
    635639    /// \param graph The digraph the algorithm runs on.
    636     /// \param arc_mixing Indicate if the arcs have to be stored in a
     640    /// \param arc_mixing Indicate if the arcs will be stored in a
    637641    /// mixed order in the internal data structure.
    638     /// In special cases, it could lead to better overall performance,
    639     /// but it is usually slower. Therefore it is disabled by default.
    640     NetworkSimplex(const GR& graph, bool arc_mixing = false) :
     642    /// In general, it leads to similar performance as using the original
     643    /// arc order, but it makes the algorithm more robust and in special
     644    /// cases, even significantly faster. Therefore, it is enabled by default.
     645    NetworkSimplex(const GR& graph, bool arc_mixing = true) :
    641646      _graph(graph), _node_id(graph), _arc_id(graph),
    642647      _arc_mixing(arc_mixing),
     
    731736    ///
    732737    /// \return <tt>(*this)</tt>
     738    ///
     739    /// \sa supplyType()
    733740    template<typename SupplyMap>
    734741    NetworkSimplex& supplyMap(const SupplyMap& map) {
     
    747754    ///
    748755    /// Using this function has the same effect as using \ref supplyMap()
    749     /// with such a map in which \c k is assigned to \c s, \c -k is
     756    /// with a map in which \c k is assigned to \c s, \c -k is
    750757    /// assigned to \c t and all other nodes have zero supply value.
    751758    ///
     
    914921      _parent.resize(all_node_num);
    915922      _pred.resize(all_node_num);
    916       _forward.resize(all_node_num);
     923      _pred_dir.resize(all_node_num);
    917924      _thread.resize(all_node_num);
    918925      _rev_thread.resize(all_node_num);
     
    928935      if (_arc_mixing) {
    929936        // Store the arcs in a mixed order
    930         int k = std::max(int(std::sqrt(double(_arc_num))), 10);
     937        const int skip = std::max(_arc_num / _node_num, 3);
    931938        int i = 0, j = 0;
    932939        for (ArcIt a(_graph); a != INVALID; ++a) {
     
    934941          _source[i] = _node_id[_graph.source(a)];
    935942          _target[i] = _node_id[_graph.target(a)];
    936           if ((i += k) >= _arc_num) i = ++j;
     943          if ((i += skip) >= _arc_num) i = ++j;
    937944        }
    938945      } else {
     
    11171124          _state[e] = STATE_TREE;
    11181125          if (_supply[u] >= 0) {
    1119             _forward[u] = true;
     1126            _pred_dir[u] = DIR_UP;
    11201127            _pi[u] = 0;
    11211128            _source[e] = u;
     
    11241131            _cost[e] = 0;
    11251132          } else {
    1126             _forward[u] = false;
     1133            _pred_dir[u] = DIR_DOWN;
    11271134            _pi[u] = ART_COST;
    11281135            _source[e] = _root;
     
    11441151          _last_succ[u] = u;
    11451152          if (_supply[u] >= 0) {
    1146             _forward[u] = true;
     1153            _pred_dir[u] = DIR_UP;
    11471154            _pi[u] = 0;
    11481155            _pred[u] = e;
     
    11541161            _state[e] = STATE_TREE;
    11551162          } else {
    1156             _forward[u] = false;
     1163            _pred_dir[u] = DIR_DOWN;
    11571164            _pi[u] = ART_COST;
    11581165            _pred[u] = f;
     
    11851192          _last_succ[u] = u;
    11861193          if (_supply[u] <= 0) {
    1187             _forward[u] = false;
     1194            _pred_dir[u] = DIR_DOWN;
    11881195            _pi[u] = 0;
    11891196            _pred[u] = e;
     
    11951202            _state[e] = STATE_TREE;
    11961203          } else {
    1197             _forward[u] = true;
     1204            _pred_dir[u] = DIR_UP;
    11981205            _pi[u] = -ART_COST;
    11991206            _pred[u] = f;
     
    12381245      // Initialize first and second nodes according to the direction
    12391246      // of the cycle
     1247      int first, second;
    12401248      if (_state[in_arc] == STATE_LOWER) {
    12411249        first  = _source[in_arc];
     
    12471255      delta = _cap[in_arc];
    12481256      int result = 0;
    1249       Value d;
     1257      Value c, d;
    12501258      int e;
    12511259
    1252       // Search the cycle along the path form the first node to the root
     1260      // Search the cycle form the first node to the join node
    12531261      for (int u = first; u != join; u = _parent[u]) {
    12541262        e = _pred[u];
    1255         d = _forward[u] ?
    1256           _flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]);
     1263        d = _flow[e];
     1264        if (_pred_dir[u] == DIR_DOWN) {
     1265          c = _cap[e];
     1266          d = c >= MAX ? INF : c - d;
     1267        }
    12571268        if (d < delta) {
    12581269          delta = d;
     
    12611272        }
    12621273      }
    1263       // Search the cycle along the path form the second node to the root
     1274
     1275      // Search the cycle form the second node to the join node
    12641276      for (int u = second; u != join; u = _parent[u]) {
    12651277        e = _pred[u];
    1266         d = _forward[u] ?
    1267           (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
     1278        d = _flow[e];
     1279        if (_pred_dir[u] == DIR_UP) {
     1280          c = _cap[e];
     1281          d = c >= MAX ? INF : c - d;
     1282        }
    12681283        if (d <= delta) {
    12691284          delta = d;
     
    12901305        _flow[in_arc] += val;
    12911306        for (int u = _source[in_arc]; u != join; u = _parent[u]) {
    1292           _flow[_pred[u]] += _forward[u] ? -val : val;
     1307          _flow[_pred[u]] -= _pred_dir[u] * val;
    12931308        }
    12941309        for (int u = _target[in_arc]; u != join; u = _parent[u]) {
    1295           _flow[_pred[u]] += _forward[u] ? val : -val;
     1310          _flow[_pred[u]] += _pred_dir[u] * val;
    12961311        }
    12971312      }
     
    13081323    // Update the tree structure
    13091324    void updateTreeStructure() {
    1310       int u, w;
    13111325      int old_rev_thread = _rev_thread[u_out];
    13121326      int old_succ_num = _succ_num[u_out];
     
    13141328      v_out = _parent[u_out];
    13151329
    1316       u = _last_succ[u_in];  // the last successor of u_in
    1317       right = _thread[u];    // the node after it
    1318 
    1319       // Handle the case when old_rev_thread equals to v_in
    1320       // (it also means that join and v_out coincide)
    1321       if (old_rev_thread == v_in) {
    1322         last = _thread[_last_succ[u_out]];
     1330      // Check if u_in and u_out coincide
     1331      if (u_in == u_out) {
     1332        // Update _parent, _pred, _pred_dir
     1333        _parent[u_in] = v_in;
     1334        _pred[u_in] = in_arc;
     1335        _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
     1336
     1337        // Update _thread and _rev_thread
     1338        if (_thread[v_in] != u_out) {
     1339          int after = _thread[old_last_succ];
     1340          _thread[old_rev_thread] = after;
     1341          _rev_thread[after] = old_rev_thread;
     1342          after = _thread[v_in];
     1343          _thread[v_in] = u_out;
     1344          _rev_thread[u_out] = v_in;
     1345          _thread[old_last_succ] = after;
     1346          _rev_thread[after] = old_last_succ;
     1347        }
    13231348      } else {
    1324         last = _thread[v_in];
    1325       }
    1326 
    1327       // Update _thread and _parent along the stem nodes (i.e. the nodes
    1328       // between u_in and u_out, whose parent have to be changed)
    1329       _thread[v_in] = stem = u_in;
    1330       _dirty_revs.clear();
    1331       _dirty_revs.push_back(v_in);
    1332       par_stem = v_in;
    1333       while (stem != u_out) {
    1334         // Insert the next stem node into the thread list
    1335         new_stem = _parent[stem];
    1336         _thread[u] = new_stem;
    1337         _dirty_revs.push_back(u);
    1338 
    1339         // Remove the subtree of stem from the thread list
    1340         w = _rev_thread[stem];
    1341         _thread[w] = right;
    1342         _rev_thread[right] = w;
    1343 
    1344         // Change the parent node and shift stem nodes
    1345         _parent[stem] = par_stem;
    1346         par_stem = stem;
    1347         stem = new_stem;
    1348 
    1349         // Update u and right
    1350         u = _last_succ[stem] == _last_succ[par_stem] ?
    1351           _rev_thread[par_stem] : _last_succ[stem];
    1352         right = _thread[u];
    1353       }
    1354       _parent[u_out] = par_stem;
    1355       _thread[u] = last;
    1356       _rev_thread[last] = u;
    1357       _last_succ[u_out] = u;
    1358 
    1359       // Remove the subtree of u_out from the thread list except for
    1360       // the case when old_rev_thread equals to v_in
    1361       // (it also means that join and v_out coincide)
    1362       if (old_rev_thread != v_in) {
    1363         _thread[old_rev_thread] = right;
    1364         _rev_thread[right] = old_rev_thread;
    1365       }
    1366 
    1367       // Update _rev_thread using the new _thread values
    1368       for (int i = 0; i != int(_dirty_revs.size()); ++i) {
    1369         u = _dirty_revs[i];
    1370         _rev_thread[_thread[u]] = u;
    1371       }
    1372 
    1373       // Update _pred, _forward, _last_succ and _succ_num for the
    1374       // stem nodes from u_out to u_in
    1375       int tmp_sc = 0, tmp_ls = _last_succ[u_out];
    1376       u = u_out;
    1377       while (u != u_in) {
    1378         w = _parent[u];
    1379         _pred[u] = _pred[w];
    1380         _forward[u] = !_forward[w];
    1381         tmp_sc += _succ_num[u] - _succ_num[w];
    1382         _succ_num[u] = tmp_sc;
    1383         _last_succ[w] = tmp_ls;
    1384         u = w;
    1385       }
    1386       _pred[u_in] = in_arc;
    1387       _forward[u_in] = (u_in == _source[in_arc]);
    1388       _succ_num[u_in] = old_succ_num;
    1389 
    1390       // Set limits for updating _last_succ form v_in and v_out
    1391       // towards the root
    1392       int up_limit_in = -1;
    1393       int up_limit_out = -1;
    1394       if (_last_succ[join] == v_in) {
    1395         up_limit_out = join;
    1396       } else {
    1397         up_limit_in = join;
     1349        // Handle the case when old_rev_thread equals to v_in
     1350        // (it also means that join and v_out coincide)
     1351        int thread_continue = old_rev_thread == v_in ?
     1352          _thread[old_last_succ] : _thread[v_in];
     1353
     1354        // Update _thread and _parent along the stem nodes (i.e. the nodes
     1355        // between u_in and u_out, whose parent have to be changed)
     1356        int stem = u_in;              // the current stem node
     1357        int par_stem = v_in;          // the new parent of stem
     1358        int next_stem;                // the next stem node
     1359        int last = _last_succ[u_in];  // the last successor of stem
     1360        int before, after = _thread[last];
     1361        _thread[v_in] = u_in;
     1362        _dirty_revs.clear();
     1363        _dirty_revs.push_back(v_in);
     1364        while (stem != u_out) {
     1365          // Insert the next stem node into the thread list
     1366          next_stem = _parent[stem];
     1367          _thread[last] = next_stem;
     1368          _dirty_revs.push_back(last);
     1369
     1370          // Remove the subtree of stem from the thread list
     1371          before = _rev_thread[stem];
     1372          _thread[before] = after;
     1373          _rev_thread[after] = before;
     1374
     1375          // Change the parent node and shift stem nodes
     1376          _parent[stem] = par_stem;
     1377          par_stem = stem;
     1378          stem = next_stem;
     1379
     1380          // Update last and after
     1381          last = _last_succ[stem] == _last_succ[par_stem] ?
     1382            _rev_thread[par_stem] : _last_succ[stem];
     1383          after = _thread[last];
     1384        }
     1385        _parent[u_out] = par_stem;
     1386        _thread[last] = thread_continue;
     1387        _rev_thread[thread_continue] = last;
     1388        _last_succ[u_out] = last;
     1389
     1390        // Remove the subtree of u_out from the thread list except for
     1391        // the case when old_rev_thread equals to v_in
     1392        if (old_rev_thread != v_in) {
     1393          _thread[old_rev_thread] = after;
     1394          _rev_thread[after] = old_rev_thread;
     1395        }
     1396
     1397        // Update _rev_thread using the new _thread values
     1398        for (int i = 0; i != int(_dirty_revs.size()); ++i) {
     1399          int u = _dirty_revs[i];
     1400          _rev_thread[_thread[u]] = u;
     1401        }
     1402
     1403        // Update _pred, _pred_dir, _last_succ and _succ_num for the
     1404        // stem nodes from u_out to u_in
     1405        int tmp_sc = 0, tmp_ls = _last_succ[u_out];
     1406        for (int u = u_out, p = _parent[u]; u != u_in; u = p, p = _parent[u]) {
     1407          _pred[u] = _pred[p];
     1408          _pred_dir[u] = -_pred_dir[p];
     1409          tmp_sc += _succ_num[u] - _succ_num[p];
     1410          _succ_num[u] = tmp_sc;
     1411          _last_succ[p] = tmp_ls;
     1412        }
     1413        _pred[u_in] = in_arc;
     1414        _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
     1415        _succ_num[u_in] = old_succ_num;
    13981416      }
    13991417
    14001418      // Update _last_succ from v_in towards the root
    1401       for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
    1402            u = _parent[u]) {
    1403         _last_succ[u] = _last_succ[u_out];
    1404       }
     1419      int up_limit_out = _last_succ[join] == v_in ? join : -1;
     1420      int last_succ_out = _last_succ[u_out];
     1421      for (int u = v_in; u != -1 && _last_succ[u] == v_in; u = _parent[u]) {
     1422        _last_succ[u] = last_succ_out;
     1423      }
     1424
    14051425      // Update _last_succ from v_out towards the root
    14061426      if (join != old_rev_thread && v_in != old_rev_thread) {
    1407         for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
     1427        for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
    14081428             u = _parent[u]) {
    14091429          _last_succ[u] = old_rev_thread;
    14101430        }
    1411       } else {
    1412         for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
     1431      }
     1432      else if (last_succ_out != old_last_succ) {
     1433        for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
    14131434             u = _parent[u]) {
    1414           _last_succ[u] = _last_succ[u_out];
     1435          _last_succ[u] = last_succ_out;
    14151436        }
    14161437      }
    14171438
    14181439      // Update _succ_num from v_in to join
    1419       for (u = v_in; u != join; u = _parent[u]) {
     1440      for (int u = v_in; u != join; u = _parent[u]) {
    14201441        _succ_num[u] += old_succ_num;
    14211442      }
    14221443      // Update _succ_num from v_out to join
    1423       for (u = v_out; u != join; u = _parent[u]) {
     1444      for (int u = v_out; u != join; u = _parent[u]) {
    14241445        _succ_num[u] -= old_succ_num;
    14251446      }
    14261447    }
    14271448
    1428     // Update potentials
     1449    // Update potentials in the subtree that has been moved
    14291450    void updatePotential() {
    1430       Cost sigma = _forward[u_in] ?
    1431         _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
    1432         _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
    1433       // Update potentials in the subtree, which has been moved
     1451      Cost sigma = _pi[v_in] - _pi[u_in] -
     1452                   _pred_dir[u_in] * _cost[in_arc];
    14341453      int end = _thread[_last_succ[u_in]];
    14351454      for (int u = u_in; u != end; u = _thread[u]) {
Note: See TracChangeset for help on using the changeset viewer.