COIN-OR::LEMON - Graph Library

Changeset 910:f3bc4e9b5f3a in lemon for lemon/network_simplex.h


Ignore:
Timestamp:
02/20/10 18:39:03 (10 years ago)
Author:
Peter Kovacs <kpeter@…>
Branch:
default
Phase:
public
Message:

New heuristics for MCF algorithms (#340)
and some implementation improvements.

  • A useful heuristic is added to NetworkSimplex? to make the initial pivots faster.
  • A powerful global update heuristic is added to CostScaling? and the implementation is reworked with various improvements.
  • Better relabeling in CostScaling? to improve numerical stability and make the code faster.
  • A small improvement is made in CapacityScaling? for better delta computation.
  • Add notes to the classes about the usage of vector<char> instead of vector<bool> for efficiency reasons.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/network_simplex.h

    r878 r910  
    165165
    166166    typedef std::vector<int> IntVector;
    167     typedef std::vector<char> CharVector;
    168167    typedef std::vector<Value> ValueVector;
    169168    typedef std::vector<Cost> CostVector;
     169    typedef std::vector<char> BoolVector;
     170    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
    170171
    171172    // State constants for arcs
     
    213214    IntVector _last_succ;
    214215    IntVector _dirty_revs;
    215     CharVector _forward;
    216     CharVector _state;
     216    BoolVector _forward;
     217    BoolVector _state;
    217218    int _root;
    218219
     
    245246      const IntVector  &_target;
    246247      const CostVector &_cost;
    247       const CharVector &_state;
     248      const BoolVector &_state;
    248249      const CostVector &_pi;
    249250      int &_in_arc;
     
    266267      bool findEnteringArc() {
    267268        Cost c;
    268         for (int e = _next_arc; e < _search_arc_num; ++e) {
     269        for (int e = _next_arc; e != _search_arc_num; ++e) {
    269270          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    270271          if (c < 0) {
     
    274275          }
    275276        }
    276         for (int e = 0; e < _next_arc; ++e) {
     277        for (int e = 0; e != _next_arc; ++e) {
    277278          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    278279          if (c < 0) {
     
    297298      const IntVector  &_target;
    298299      const CostVector &_cost;
    299       const CharVector &_state;
     300      const BoolVector &_state;
    300301      const CostVector &_pi;
    301302      int &_in_arc;
     
    314315      bool findEnteringArc() {
    315316        Cost c, min = 0;
    316         for (int e = 0; e < _search_arc_num; ++e) {
     317        for (int e = 0; e != _search_arc_num; ++e) {
    317318          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    318319          if (c < min) {
     
    336337      const IntVector  &_target;
    337338      const CostVector &_cost;
    338       const CharVector &_state;
     339      const BoolVector &_state;
    339340      const CostVector &_pi;
    340341      int &_in_arc;
     
    355356      {
    356357        // The main parameters of the pivot rule
    357         const double BLOCK_SIZE_FACTOR = 0.5;
     358        const double BLOCK_SIZE_FACTOR = 1.0;
    358359        const int MIN_BLOCK_SIZE = 10;
    359360
     
    368369        int cnt = _block_size;
    369370        int e;
    370         for (e = _next_arc; e < _search_arc_num; ++e) {
     371        for (e = _next_arc; e != _search_arc_num; ++e) {
    371372          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    372373          if (c < min) {
     
    379380          }
    380381        }
    381         for (e = 0; e < _next_arc; ++e) {
     382        for (e = 0; e != _next_arc; ++e) {
    382383          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    383384          if (c < min) {
     
    409410      const IntVector  &_target;
    410411      const CostVector &_cost;
    411       const CharVector &_state;
     412      const BoolVector &_state;
    412413      const CostVector &_pi;
    413414      int &_in_arc;
     
    470471        min = 0;
    471472        _curr_length = 0;
    472         for (e = _next_arc; e < _search_arc_num; ++e) {
     473        for (e = _next_arc; e != _search_arc_num; ++e) {
    473474          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    474475          if (c < 0) {
     
    481482          }
    482483        }
    483         for (e = 0; e < _next_arc; ++e) {
     484        for (e = 0; e != _next_arc; ++e) {
    484485          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    485486          if (c < 0) {
     
    512513      const IntVector  &_target;
    513514      const CostVector &_cost;
    514       const CharVector &_state;
     515      const BoolVector &_state;
    515516      const CostVector &_pi;
    516517      int &_in_arc;
     
    565566        // Check the current candidate list
    566567        int e;
    567         for (int i = 0; i < _curr_length; ++i) {
     568        for (int i = 0; i != _curr_length; ++i) {
    568569          e = _candidates[i];
    569570          _cand_cost[e] = _state[e] *
     
    578579        int limit = _head_length;
    579580
    580         for (e = _next_arc; e < _search_arc_num; ++e) {
     581        for (e = _next_arc; e != _search_arc_num; ++e) {
    581582          _cand_cost[e] = _state[e] *
    582583            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     
    590591          }
    591592        }
    592         for (e = 0; e < _next_arc; ++e) {
     593        for (e = 0; e != _next_arc; ++e) {
    593594          _cand_cost[e] = _state[e] *
    594595            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     
    13291330
    13301331      // Update _rev_thread using the new _thread values
    1331       for (int i = 0; i < int(_dirty_revs.size()); ++i) {
     1332      for (int i = 0; i != int(_dirty_revs.size()); ++i) {
    13321333        u = _dirty_revs[i];
    13331334        _rev_thread[_thread[u]] = u;
     
    13991400        _pi[u] += sigma;
    14001401      }
     1402    }
     1403
     1404    // Heuristic initial pivots
     1405    bool initialPivots() {
     1406      Value curr, total = 0;
     1407      std::vector<Node> supply_nodes, demand_nodes;
     1408      for (NodeIt u(_graph); u != INVALID; ++u) {
     1409        curr = _supply[_node_id[u]];
     1410        if (curr > 0) {
     1411          total += curr;
     1412          supply_nodes.push_back(u);
     1413        }
     1414        else if (curr < 0) {
     1415          demand_nodes.push_back(u);
     1416        }
     1417      }
     1418      if (_sum_supply > 0) total -= _sum_supply;
     1419      if (total <= 0) return true;
     1420
     1421      IntVector arc_vector;
     1422      if (_sum_supply >= 0) {
     1423        if (supply_nodes.size() == 1 && demand_nodes.size() == 1) {
     1424          // Perform a reverse graph search from the sink to the source
     1425          typename GR::template NodeMap<bool> reached(_graph, false);
     1426          Node s = supply_nodes[0], t = demand_nodes[0];
     1427          std::vector<Node> stack;
     1428          reached[t] = true;
     1429          stack.push_back(t);
     1430          while (!stack.empty()) {
     1431            Node u, v = stack.back();
     1432            stack.pop_back();
     1433            if (v == s) break;
     1434            for (InArcIt a(_graph, v); a != INVALID; ++a) {
     1435              if (reached[u = _graph.source(a)]) continue;
     1436              int j = _arc_id[a];
     1437              if (_cap[j] >= total) {
     1438                arc_vector.push_back(j);
     1439                reached[u] = true;
     1440                stack.push_back(u);
     1441              }
     1442            }
     1443          }
     1444        } else {
     1445          // Find the min. cost incomming arc for each demand node
     1446          for (int i = 0; i != int(demand_nodes.size()); ++i) {
     1447            Node v = demand_nodes[i];
     1448            Cost c, min_cost = std::numeric_limits<Cost>::max();
     1449            Arc min_arc = INVALID;
     1450            for (InArcIt a(_graph, v); a != INVALID; ++a) {
     1451              c = _cost[_arc_id[a]];
     1452              if (c < min_cost) {
     1453                min_cost = c;
     1454                min_arc = a;
     1455              }
     1456            }
     1457            if (min_arc != INVALID) {
     1458              arc_vector.push_back(_arc_id[min_arc]);
     1459            }
     1460          }
     1461        }
     1462      } else {
     1463        // Find the min. cost outgoing arc for each supply node
     1464        for (int i = 0; i != int(supply_nodes.size()); ++i) {
     1465          Node u = supply_nodes[i];
     1466          Cost c, min_cost = std::numeric_limits<Cost>::max();
     1467          Arc min_arc = INVALID;
     1468          for (OutArcIt a(_graph, u); a != INVALID; ++a) {
     1469            c = _cost[_arc_id[a]];
     1470            if (c < min_cost) {
     1471              min_cost = c;
     1472              min_arc = a;
     1473            }
     1474          }
     1475          if (min_arc != INVALID) {
     1476            arc_vector.push_back(_arc_id[min_arc]);
     1477          }
     1478        }
     1479      }
     1480
     1481      // Perform heuristic initial pivots
     1482      for (int i = 0; i != int(arc_vector.size()); ++i) {
     1483        in_arc = arc_vector[i];
     1484        if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] -
     1485            _pi[_target[in_arc]]) >= 0) continue;
     1486        findJoinNode();
     1487        bool change = findLeavingArc();
     1488        if (delta >= MAX) return false;
     1489        changeFlow(change);
     1490        if (change) {
     1491          updateTreeStructure();
     1492          updatePotential();
     1493        }
     1494      }
     1495      return true;
    14011496    }
    14021497
     
    14231518      PivotRuleImpl pivot(*this);
    14241519
     1520      // Perform heuristic initial pivots
     1521      if (!initialPivots()) return UNBOUNDED;
     1522
    14251523      // Execute the Network Simplex algorithm
    14261524      while (pivot.findEnteringArc()) {
Note: See TracChangeset for help on using the changeset viewer.