lemon/network_simplex.h
changeset 857 abb95d48e89e
parent 830 75c97c3786d6
parent 839 f3bc4e9b5f3a
child 862 b6f76c95992e
equal deleted inserted replaced
25:29f2d3636085 27:b74e7ff13cc7
   162   private:
   162   private:
   163 
   163 
   164     TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   164     TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   165 
   165 
   166     typedef std::vector<int> IntVector;
   166     typedef std::vector<int> IntVector;
   167     typedef std::vector<char> CharVector;
       
   168     typedef std::vector<Value> ValueVector;
   167     typedef std::vector<Value> ValueVector;
   169     typedef std::vector<Cost> CostVector;
   168     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 
   171 
   171     // State constants for arcs
   172     // State constants for arcs
   172     enum ArcStateEnum {
   173     enum ArcStateEnum {
   173       STATE_UPPER = -1,
   174       STATE_UPPER = -1,
   174       STATE_TREE  =  0,
   175       STATE_TREE  =  0,
   211     IntVector _thread;
   212     IntVector _thread;
   212     IntVector _rev_thread;
   213     IntVector _rev_thread;
   213     IntVector _succ_num;
   214     IntVector _succ_num;
   214     IntVector _last_succ;
   215     IntVector _last_succ;
   215     IntVector _dirty_revs;
   216     IntVector _dirty_revs;
   216     CharVector _forward;
   217     BoolVector _forward;
   217     CharVector _state;
   218     BoolVector _state;
   218     int _root;
   219     int _root;
   219 
   220 
   220     // Temporary data used in the current pivot iteration
   221     // Temporary data used in the current pivot iteration
   221     int in_arc, join, u_in, v_in, u_out, v_out;
   222     int in_arc, join, u_in, v_in, u_out, v_out;
   222     int first, second, right, last;
   223     int first, second, right, last;
   243 
   244 
   244       // References to the NetworkSimplex class
   245       // References to the NetworkSimplex class
   245       const IntVector  &_source;
   246       const IntVector  &_source;
   246       const IntVector  &_target;
   247       const IntVector  &_target;
   247       const CostVector &_cost;
   248       const CostVector &_cost;
   248       const CharVector &_state;
   249       const BoolVector &_state;
   249       const CostVector &_pi;
   250       const CostVector &_pi;
   250       int &_in_arc;
   251       int &_in_arc;
   251       int _search_arc_num;
   252       int _search_arc_num;
   252 
   253 
   253       // Pivot rule data
   254       // Pivot rule data
   264       {}
   265       {}
   265 
   266 
   266       // Find next entering arc
   267       // Find next entering arc
   267       bool findEnteringArc() {
   268       bool findEnteringArc() {
   268         Cost c;
   269         Cost c;
   269         for (int e = _next_arc; e < _search_arc_num; ++e) {
   270         for (int e = _next_arc; e != _search_arc_num; ++e) {
   270           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   271           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   271           if (c < 0) {
   272           if (c < 0) {
   272             _in_arc = e;
   273             _in_arc = e;
   273             _next_arc = e + 1;
   274             _next_arc = e + 1;
   274             return true;
   275             return true;
   275           }
   276           }
   276         }
   277         }
   277         for (int e = 0; e < _next_arc; ++e) {
   278         for (int e = 0; e != _next_arc; ++e) {
   278           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   279           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   279           if (c < 0) {
   280           if (c < 0) {
   280             _in_arc = e;
   281             _in_arc = e;
   281             _next_arc = e + 1;
   282             _next_arc = e + 1;
   282             return true;
   283             return true;
   295 
   296 
   296       // References to the NetworkSimplex class
   297       // References to the NetworkSimplex class
   297       const IntVector  &_source;
   298       const IntVector  &_source;
   298       const IntVector  &_target;
   299       const IntVector  &_target;
   299       const CostVector &_cost;
   300       const CostVector &_cost;
   300       const CharVector &_state;
   301       const BoolVector &_state;
   301       const CostVector &_pi;
   302       const CostVector &_pi;
   302       int &_in_arc;
   303       int &_in_arc;
   303       int _search_arc_num;
   304       int _search_arc_num;
   304 
   305 
   305     public:
   306     public:
   312       {}
   313       {}
   313 
   314 
   314       // Find next entering arc
   315       // Find next entering arc
   315       bool findEnteringArc() {
   316       bool findEnteringArc() {
   316         Cost c, min = 0;
   317         Cost c, min = 0;
   317         for (int e = 0; e < _search_arc_num; ++e) {
   318         for (int e = 0; e != _search_arc_num; ++e) {
   318           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   319           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   319           if (c < min) {
   320           if (c < min) {
   320             min = c;
   321             min = c;
   321             _in_arc = e;
   322             _in_arc = e;
   322           }
   323           }
   334 
   335 
   335       // References to the NetworkSimplex class
   336       // References to the NetworkSimplex class
   336       const IntVector  &_source;
   337       const IntVector  &_source;
   337       const IntVector  &_target;
   338       const IntVector  &_target;
   338       const CostVector &_cost;
   339       const CostVector &_cost;
   339       const CharVector &_state;
   340       const BoolVector &_state;
   340       const CostVector &_pi;
   341       const CostVector &_pi;
   341       int &_in_arc;
   342       int &_in_arc;
   342       int _search_arc_num;
   343       int _search_arc_num;
   343 
   344 
   344       // Pivot rule data
   345       // Pivot rule data
   353         _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   354         _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   354         _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
   355         _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
   355         _next_arc(0)
   356         _next_arc(0)
   356       {
   357       {
   357         // The main parameters of the pivot rule
   358         // The main parameters of the pivot rule
   358         const double BLOCK_SIZE_FACTOR = 0.5;
   359         const double BLOCK_SIZE_FACTOR = 1.0;
   359         const int MIN_BLOCK_SIZE = 10;
   360         const int MIN_BLOCK_SIZE = 10;
   360 
   361 
   361         _block_size = std::max( int(BLOCK_SIZE_FACTOR *
   362         _block_size = std::max( int(BLOCK_SIZE_FACTOR *
   362                                     std::sqrt(double(_search_arc_num))),
   363                                     std::sqrt(double(_search_arc_num))),
   363                                 MIN_BLOCK_SIZE );
   364                                 MIN_BLOCK_SIZE );
   366       // Find next entering arc
   367       // Find next entering arc
   367       bool findEnteringArc() {
   368       bool findEnteringArc() {
   368         Cost c, min = 0;
   369         Cost c, min = 0;
   369         int cnt = _block_size;
   370         int cnt = _block_size;
   370         int e;
   371         int e;
   371         for (e = _next_arc; e < _search_arc_num; ++e) {
   372         for (e = _next_arc; e != _search_arc_num; ++e) {
   372           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   373           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   373           if (c < min) {
   374           if (c < min) {
   374             min = c;
   375             min = c;
   375             _in_arc = e;
   376             _in_arc = e;
   376           }
   377           }
   377           if (--cnt == 0) {
   378           if (--cnt == 0) {
   378             if (min < 0) goto search_end;
   379             if (min < 0) goto search_end;
   379             cnt = _block_size;
   380             cnt = _block_size;
   380           }
   381           }
   381         }
   382         }
   382         for (e = 0; e < _next_arc; ++e) {
   383         for (e = 0; e != _next_arc; ++e) {
   383           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   384           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   384           if (c < min) {
   385           if (c < min) {
   385             min = c;
   386             min = c;
   386             _in_arc = e;
   387             _in_arc = e;
   387           }
   388           }
   407 
   408 
   408       // References to the NetworkSimplex class
   409       // References to the NetworkSimplex class
   409       const IntVector  &_source;
   410       const IntVector  &_source;
   410       const IntVector  &_target;
   411       const IntVector  &_target;
   411       const CostVector &_cost;
   412       const CostVector &_cost;
   412       const CharVector &_state;
   413       const BoolVector &_state;
   413       const CostVector &_pi;
   414       const CostVector &_pi;
   414       int &_in_arc;
   415       int &_in_arc;
   415       int _search_arc_num;
   416       int _search_arc_num;
   416 
   417 
   417       // Pivot rule data
   418       // Pivot rule data
   468         }
   469         }
   469 
   470 
   470         // Major iteration: build a new candidate list
   471         // Major iteration: build a new candidate list
   471         min = 0;
   472         min = 0;
   472         _curr_length = 0;
   473         _curr_length = 0;
   473         for (e = _next_arc; e < _search_arc_num; ++e) {
   474         for (e = _next_arc; e != _search_arc_num; ++e) {
   474           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   475           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   475           if (c < 0) {
   476           if (c < 0) {
   476             _candidates[_curr_length++] = e;
   477             _candidates[_curr_length++] = e;
   477             if (c < min) {
   478             if (c < min) {
   478               min = c;
   479               min = c;
   479               _in_arc = e;
   480               _in_arc = e;
   480             }
   481             }
   481             if (_curr_length == _list_length) goto search_end;
   482             if (_curr_length == _list_length) goto search_end;
   482           }
   483           }
   483         }
   484         }
   484         for (e = 0; e < _next_arc; ++e) {
   485         for (e = 0; e != _next_arc; ++e) {
   485           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   486           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   486           if (c < 0) {
   487           if (c < 0) {
   487             _candidates[_curr_length++] = e;
   488             _candidates[_curr_length++] = e;
   488             if (c < min) {
   489             if (c < min) {
   489               min = c;
   490               min = c;
   510 
   511 
   511       // References to the NetworkSimplex class
   512       // References to the NetworkSimplex class
   512       const IntVector  &_source;
   513       const IntVector  &_source;
   513       const IntVector  &_target;
   514       const IntVector  &_target;
   514       const CostVector &_cost;
   515       const CostVector &_cost;
   515       const CharVector &_state;
   516       const BoolVector &_state;
   516       const CostVector &_pi;
   517       const CostVector &_pi;
   517       int &_in_arc;
   518       int &_in_arc;
   518       int _search_arc_num;
   519       int _search_arc_num;
   519 
   520 
   520       // Pivot rule data
   521       // Pivot rule data
   563 
   564 
   564       // Find next entering arc
   565       // Find next entering arc
   565       bool findEnteringArc() {
   566       bool findEnteringArc() {
   566         // Check the current candidate list
   567         // Check the current candidate list
   567         int e;
   568         int e;
   568         for (int i = 0; i < _curr_length; ++i) {
   569         for (int i = 0; i != _curr_length; ++i) {
   569           e = _candidates[i];
   570           e = _candidates[i];
   570           _cand_cost[e] = _state[e] *
   571           _cand_cost[e] = _state[e] *
   571             (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   572             (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   572           if (_cand_cost[e] >= 0) {
   573           if (_cand_cost[e] >= 0) {
   573             _candidates[i--] = _candidates[--_curr_length];
   574             _candidates[i--] = _candidates[--_curr_length];
   576 
   577 
   577         // Extend the list
   578         // Extend the list
   578         int cnt = _block_size;
   579         int cnt = _block_size;
   579         int limit = _head_length;
   580         int limit = _head_length;
   580 
   581 
   581         for (e = _next_arc; e < _search_arc_num; ++e) {
   582         for (e = _next_arc; e != _search_arc_num; ++e) {
   582           _cand_cost[e] = _state[e] *
   583           _cand_cost[e] = _state[e] *
   583             (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   584             (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   584           if (_cand_cost[e] < 0) {
   585           if (_cand_cost[e] < 0) {
   585             _candidates[_curr_length++] = e;
   586             _candidates[_curr_length++] = e;
   586           }
   587           }
   588             if (_curr_length > limit) goto search_end;
   589             if (_curr_length > limit) goto search_end;
   589             limit = 0;
   590             limit = 0;
   590             cnt = _block_size;
   591             cnt = _block_size;
   591           }
   592           }
   592         }
   593         }
   593         for (e = 0; e < _next_arc; ++e) {
   594         for (e = 0; e != _next_arc; ++e) {
   594           _cand_cost[e] = _state[e] *
   595           _cand_cost[e] = _state[e] *
   595             (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   596             (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   596           if (_cand_cost[e] < 0) {
   597           if (_cand_cost[e] < 0) {
   597             _candidates[_curr_length++] = e;
   598             _candidates[_curr_length++] = e;
   598           }
   599           }
  1358         _thread[old_rev_thread] = right;
  1359         _thread[old_rev_thread] = right;
  1359         _rev_thread[right] = old_rev_thread;
  1360         _rev_thread[right] = old_rev_thread;
  1360       }
  1361       }
  1361 
  1362 
  1362       // Update _rev_thread using the new _thread values
  1363       // Update _rev_thread using the new _thread values
  1363       for (int i = 0; i < int(_dirty_revs.size()); ++i) {
  1364       for (int i = 0; i != int(_dirty_revs.size()); ++i) {
  1364         u = _dirty_revs[i];
  1365         u = _dirty_revs[i];
  1365         _rev_thread[_thread[u]] = u;
  1366         _rev_thread[_thread[u]] = u;
  1366       }
  1367       }
  1367 
  1368 
  1368       // Update _pred, _forward, _last_succ and _succ_num for the
  1369       // Update _pred, _forward, _last_succ and _succ_num for the
  1428       // Update potentials in the subtree, which has been moved
  1429       // Update potentials in the subtree, which has been moved
  1429       int end = _thread[_last_succ[u_in]];
  1430       int end = _thread[_last_succ[u_in]];
  1430       for (int u = u_in; u != end; u = _thread[u]) {
  1431       for (int u = u_in; u != end; u = _thread[u]) {
  1431         _pi[u] += sigma;
  1432         _pi[u] += sigma;
  1432       }
  1433       }
       
  1434     }
       
  1435 
       
  1436     // Heuristic initial pivots
       
  1437     bool initialPivots() {
       
  1438       Value curr, total = 0;
       
  1439       std::vector<Node> supply_nodes, demand_nodes;
       
  1440       for (NodeIt u(_graph); u != INVALID; ++u) {
       
  1441         curr = _supply[_node_id[u]];
       
  1442         if (curr > 0) {
       
  1443           total += curr;
       
  1444           supply_nodes.push_back(u);
       
  1445         }
       
  1446         else if (curr < 0) {
       
  1447           demand_nodes.push_back(u);
       
  1448         }
       
  1449       }
       
  1450       if (_sum_supply > 0) total -= _sum_supply;
       
  1451       if (total <= 0) return true;
       
  1452 
       
  1453       IntVector arc_vector;
       
  1454       if (_sum_supply >= 0) {
       
  1455         if (supply_nodes.size() == 1 && demand_nodes.size() == 1) {
       
  1456           // Perform a reverse graph search from the sink to the source
       
  1457           typename GR::template NodeMap<bool> reached(_graph, false);
       
  1458           Node s = supply_nodes[0], t = demand_nodes[0];
       
  1459           std::vector<Node> stack;
       
  1460           reached[t] = true;
       
  1461           stack.push_back(t);
       
  1462           while (!stack.empty()) {
       
  1463             Node u, v = stack.back();
       
  1464             stack.pop_back();
       
  1465             if (v == s) break;
       
  1466             for (InArcIt a(_graph, v); a != INVALID; ++a) {
       
  1467               if (reached[u = _graph.source(a)]) continue;
       
  1468               int j = _arc_id[a];
       
  1469               if (_cap[j] >= total) {
       
  1470                 arc_vector.push_back(j);
       
  1471                 reached[u] = true;
       
  1472                 stack.push_back(u);
       
  1473               }
       
  1474             }
       
  1475           }
       
  1476         } else {
       
  1477           // Find the min. cost incomming arc for each demand node
       
  1478           for (int i = 0; i != int(demand_nodes.size()); ++i) {
       
  1479             Node v = demand_nodes[i];
       
  1480             Cost c, min_cost = std::numeric_limits<Cost>::max();
       
  1481             Arc min_arc = INVALID;
       
  1482             for (InArcIt a(_graph, v); a != INVALID; ++a) {
       
  1483               c = _cost[_arc_id[a]];
       
  1484               if (c < min_cost) {
       
  1485                 min_cost = c;
       
  1486                 min_arc = a;
       
  1487               }
       
  1488             }
       
  1489             if (min_arc != INVALID) {
       
  1490               arc_vector.push_back(_arc_id[min_arc]);
       
  1491             }
       
  1492           }
       
  1493         }
       
  1494       } else {
       
  1495         // Find the min. cost outgoing arc for each supply node
       
  1496         for (int i = 0; i != int(supply_nodes.size()); ++i) {
       
  1497           Node u = supply_nodes[i];
       
  1498           Cost c, min_cost = std::numeric_limits<Cost>::max();
       
  1499           Arc min_arc = INVALID;
       
  1500           for (OutArcIt a(_graph, u); a != INVALID; ++a) {
       
  1501             c = _cost[_arc_id[a]];
       
  1502             if (c < min_cost) {
       
  1503               min_cost = c;
       
  1504               min_arc = a;
       
  1505             }
       
  1506           }
       
  1507           if (min_arc != INVALID) {
       
  1508             arc_vector.push_back(_arc_id[min_arc]);
       
  1509           }
       
  1510         }
       
  1511       }
       
  1512 
       
  1513       // Perform heuristic initial pivots
       
  1514       for (int i = 0; i != int(arc_vector.size()); ++i) {
       
  1515         in_arc = arc_vector[i];
       
  1516         if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] -
       
  1517             _pi[_target[in_arc]]) >= 0) continue;
       
  1518         findJoinNode();
       
  1519         bool change = findLeavingArc();
       
  1520         if (delta >= MAX) return false;
       
  1521         changeFlow(change);
       
  1522         if (change) {
       
  1523           updateTreeStructure();
       
  1524           updatePotential();
       
  1525         }
       
  1526       }
       
  1527       return true;
  1433     }
  1528     }
  1434 
  1529 
  1435     // Execute the algorithm
  1530     // Execute the algorithm
  1436     ProblemType start(PivotRule pivot_rule) {
  1531     ProblemType start(PivotRule pivot_rule) {
  1437       // Select the pivot rule implementation
  1532       // Select the pivot rule implementation
  1452 
  1547 
  1453     template <typename PivotRuleImpl>
  1548     template <typename PivotRuleImpl>
  1454     ProblemType start() {
  1549     ProblemType start() {
  1455       PivotRuleImpl pivot(*this);
  1550       PivotRuleImpl pivot(*this);
  1456 
  1551 
       
  1552       // Perform heuristic initial pivots
       
  1553       if (!initialPivots()) return UNBOUNDED;
       
  1554 
  1457       // Execute the Network Simplex algorithm
  1555       // Execute the Network Simplex algorithm
  1458       while (pivot.findEnteringArc()) {
  1556       while (pivot.findEnteringArc()) {
  1459         findJoinNode();
  1557         findJoinNode();
  1460         bool change = findLeavingArc();
  1558         bool change = findLeavingArc();
  1461         if (delta >= MAX) return UNBOUNDED;
  1559         if (delta >= MAX) return UNBOUNDED;