COIN-OR::LEMON - Graph Library

Changeset 2575:e866e288cba6 in lemon-0.x


Ignore:
Timestamp:
02/18/08 04:32:06 (11 years ago)
Author:
Peter Kovacs
Branch:
default
Phase:
public
Convert:
svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3457
Message:

Major improvements in NetworkSimplex?.

Main changes:

  • Use -potenital[] instead of potential[] to conform to the usual

terminology.

  • Use function parameter instead of #define commands to select pivot rule.
  • Use much faster implementation for the candidate list pivot rule.

It is about 5-20 times faster now.

  • Add a new pivot rule called "Limited Search" that is a modified

version of "Block Search". It is about 25 percent faster on rather
sparse graphs.

  • By default "Limited Search" is used for sparse graphs and

"Block Search" is used otherwise. This combined method is the most
efficient on every input class.

  • Change the name of private members to start with "_".
  • Change the name of function parameters not to start with "_".
  • Remove unnecessary documentation for private members.
  • Many doc improvements.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/network_simplex.h

    r2569 r2575  
    2323///
    2424/// \file
    25 /// \brief The network simplex algorithm for finding a minimum cost flow.
    26 
     25/// \brief Network simplex algorithm for finding a minimum cost flow.
     26
     27#include <vector>
    2728#include <limits>
     29
    2830#include <lemon/graph_adaptor.h>
    2931#include <lemon/graph_utils.h>
    3032#include <lemon/smart_graph.h>
    31 
    32 /// \brief The pivot rule used in the algorithm.
    33 //#define FIRST_ELIGIBLE_PIVOT
    34 //#define BEST_ELIGIBLE_PIVOT
    35 #define EDGE_BLOCK_PIVOT
    36 //#define CANDIDATE_LIST_PIVOT
    37 //#define SORTED_LIST_PIVOT
    38 
    39 //#define _DEBUG_ITER_
    40 
    41 
    42 // State constant for edges at their lower bounds.
    43 #define LOWER   1
    44 // State constant for edges in the spanning tree.
    45 #define TREE    0
    46 // State constant for edges at their upper bounds.
    47 #define UPPER   -1
    48 
    49 #ifdef EDGE_BLOCK_PIVOT
    50   #include <lemon/math.h>
    51   #define MIN_BLOCK_SIZE        10
    52 #endif
    53 
    54 #ifdef CANDIDATE_LIST_PIVOT
    55   #include <vector>
    56   #define LIST_LENGTH_DIV       50
    57   #define MINOR_LIMIT_DIV       200
    58 #endif
    59 
    60 #ifdef SORTED_LIST_PIVOT
    61   #include <vector>
    62   #include <algorithm>
    63   #define LIST_LENGTH_DIV       100
    64   #define LOWER_DIV             4
    65 #endif
     33#include <lemon/math.h>
    6634
    6735namespace lemon {
     
    7644  /// finding a minimum cost flow.
    7745  ///
    78   /// \param Graph The directed graph type the algorithm runs on.
    79   /// \param LowerMap The type of the lower bound map.
    80   /// \param CapacityMap The type of the capacity (upper bound) map.
    81   /// \param CostMap The type of the cost (length) map.
    82   /// \param SupplyMap The type of the supply map.
     46  /// \tparam Graph The directed graph type the algorithm runs on.
     47  /// \tparam LowerMap The type of the lower bound map.
     48  /// \tparam CapacityMap The type of the capacity (upper bound) map.
     49  /// \tparam CostMap The type of the cost (length) map.
     50  /// \tparam SupplyMap The type of the supply map.
    8351  ///
    8452  /// \warning
    85   /// - Edge capacities and costs should be non-negative integers.
    86   ///   However \c CostMap::Value should be signed type.
    87   /// - Supply values should be signed integers.
    88   /// - \c LowerMap::Value must be convertible to
    89   ///   \c CapacityMap::Value and \c CapacityMap::Value must be
    90   ///   convertible to \c SupplyMap::Value.
     53  /// - Edge capacities and costs should be \e non-negative \e integers.
     54  /// - Supply values should be \e signed \e integers.
     55  /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.
     56  /// - \c CapacityMap::Value and \c SupplyMap::Value must be
     57  ///   convertible to each other.
     58  /// - All value types must be convertible to \c CostMap::Value, which
     59  ///   must be signed type.
     60  ///
     61  /// \note \ref NetworkSimplex provides six different pivot rule
     62  /// implementations that significantly affect the efficiency of the
     63  /// algorithm.
     64  /// By default a combined pivot rule is used, which is the fastest
     65  /// implementation according to our benchmark tests.
     66  /// Another pivot rule can be selected using \ref run() function
     67  /// with the proper parameter.
    9168  ///
    9269  /// \author Peter Kovacs
     
    9471  template < typename Graph,
    9572             typename LowerMap = typename Graph::template EdgeMap<int>,
    96              typename CapacityMap = LowerMap,
     73             typename CapacityMap = typename Graph::template EdgeMap<int>,
    9774             typename CostMap = typename Graph::template EdgeMap<int>,
    98              typename SupplyMap = typename Graph::template NodeMap
    99                                   <typename CapacityMap::Value> >
     75             typename SupplyMap = typename Graph::template NodeMap<int> >
    10076  class NetworkSimplex
    10177  {
    102     typedef typename LowerMap::Value Lower;
    10378    typedef typename CapacityMap::Value Capacity;
    10479    typedef typename CostMap::Value Cost;
     
    10883    GRAPH_TYPEDEFS(typename SGraph);
    10984
    110     typedef typename SGraph::template EdgeMap<Lower> SLowerMap;
    11185    typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap;
    11286    typedef typename SGraph::template EdgeMap<Cost> SCostMap;
     
    130104    typedef typename Graph::template NodeMap<Cost> PotentialMap;
    131105
    132   protected:
    133 
     106  public:
     107
     108    /// Enum type to select the pivot rule used by \ref run().
     109    enum PivotRuleEnum {
     110      FIRST_ELIGIBLE_PIVOT,
     111      BEST_ELIGIBLE_PIVOT,
     112      BLOCK_SEARCH_PIVOT,
     113      LIMITED_SEARCH_PIVOT,
     114      CANDIDATE_LIST_PIVOT,
     115      COMBINED_PIVOT
     116    };
     117
     118  private:
     119
     120    /// \brief Map adaptor class for handling reduced edge costs.
     121    ///
    134122    /// Map adaptor class for handling reduced edge costs.
    135123    class ReducedCostMap : public MapBase<Edge, Cost>
     
    137125    private:
    138126
    139       const SGraph &gr;
    140       const SCostMap &cost_map;
    141       const SPotentialMap &pot_map;
     127      const SGraph &_gr;
     128      const SCostMap &_cost_map;
     129      const SPotentialMap &_pot_map;
    142130
    143131    public:
    144132
    145       ReducedCostMap( const SGraph &_gr,
    146                       const SCostMap &_cm,
    147                       const SPotentialMap &_pm ) :
    148         gr(_gr), cost_map(_cm), pot_map(_pm) {}
    149 
     133      ///\e
     134      ReducedCostMap( const SGraph &gr,
     135                      const SCostMap &cost_map,
     136                      const SPotentialMap &pot_map ) :
     137        _gr(gr), _cost_map(cost_map), _pot_map(pm) {}
     138
     139      ///\e
    150140      Cost operator[](const Edge &e) const {
    151         return cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];
     141        return _cost_map[e] + _pot_map[_gr.source(e)]
     142                           - _pot_map[_gr.target(e)];
    152143      }
    153144
    154145    }; //class ReducedCostMap
    155146
    156   protected:
    157 
    158     /// The directed graph the algorithm runs on.
    159     SGraph graph;
    160     /// The original graph.
    161     const Graph &graph_ref;
    162     /// The original lower bound map.
    163     const LowerMap *lower;
    164     /// The capacity map.
    165     SCapacityMap capacity;
    166     /// The cost map.
    167     SCostMap cost;
    168     /// The supply map.
    169     SSupplyMap supply;
    170     /// The reduced cost map.
    171     ReducedCostMap red_cost;
    172     bool valid_supply;
    173 
    174     /// The edge map of the current flow.
    175     SCapacityMap flow;
    176     /// The potential node map.
    177     SPotentialMap potential;
    178     FlowMap flow_result;
    179     PotentialMap potential_result;
    180 
    181     /// Node reference for the original graph.
    182     NodeRefMap node_ref;
    183     /// Edge reference for the original graph.
    184     EdgeRefMap edge_ref;
    185 
    186     /// The \c depth node map of the spanning tree structure.
    187     IntNodeMap depth;
    188     /// The \c parent node map of the spanning tree structure.
    189     NodeNodeMap parent;
    190     /// The \c pred_edge node map of the spanning tree structure.
    191     EdgeNodeMap pred_edge;
    192     /// The \c thread node map of the spanning tree structure.
    193     NodeNodeMap thread;
    194     /// The \c forward node map of the spanning tree structure.
    195     BoolNodeMap forward;
    196     /// The state edge map.
    197     IntEdgeMap state;
    198     /// The root node of the starting spanning tree.
    199     Node root;
     147  private:
     148
     149    /// \brief Implementation of the "First Eligible" pivot rule for the
     150    /// \ref NetworkSimplex "network simplex" algorithm.
     151    ///
     152    /// This class implements the "First Eligible" pivot rule
     153    /// for the \ref NetworkSimplex "network simplex" algorithm.
     154    class FirstEligiblePivotRule
     155    {
     156    private:
     157
     158      NetworkSimplex &_ns;
     159      EdgeIt _next_edge;
     160
     161    public:
     162
     163      /// Constructor.
     164      FirstEligiblePivotRule(NetworkSimplex &ns) :
     165        _ns(ns), _next_edge(ns._graph) {}
     166
     167      /// Finds the next entering edge.
     168      bool findEnteringEdge() {
     169        for (EdgeIt e = _next_edge; e != INVALID; ++e) {
     170          if (_ns._state[e] * _ns._red_cost[e] < 0) {
     171            _ns._in_edge = e;
     172            _next_edge = ++e;
     173            return true;
     174          }
     175        }
     176        for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
     177          if (_ns._state[e] * _ns._red_cost[e] < 0) {
     178            _ns._in_edge = e;
     179            _next_edge = ++e;
     180            return true;
     181          }
     182        }
     183        return false;
     184      }
     185    }; //class FirstEligiblePivotRule
     186
     187    /// \brief Implementation of the "Best Eligible" pivot rule for the
     188    /// \ref NetworkSimplex "network simplex" algorithm.
     189    ///
     190    /// This class implements the "Best Eligible" pivot rule
     191    /// for the \ref NetworkSimplex "network simplex" algorithm.
     192    class BestEligiblePivotRule
     193    {
     194    private:
     195
     196      NetworkSimplex &_ns;
     197
     198    public:
     199
     200      /// Constructor.
     201      BestEligiblePivotRule(NetworkSimplex &ns) : _ns(ns) {}
     202
     203      /// Finds the next entering edge.
     204      bool findEnteringEdge() {
     205        Cost min = 0;
     206        for (EdgeIt e(_ns._graph); e != INVALID; ++e) {
     207          if (_ns._state[e] * _ns._red_cost[e] < min) {
     208            min = _ns._state[e] * _ns._red_cost[e];
     209            _ns._in_edge = e;
     210          }
     211        }
     212        return min < 0;
     213      }
     214    }; //class BestEligiblePivotRule
     215
     216    /// \brief Implementation of the "Block Search" pivot rule for the
     217    /// \ref NetworkSimplex "network simplex" algorithm.
     218    ///
     219    /// This class implements the "Block Search" pivot rule
     220    /// for the \ref NetworkSimplex "network simplex" algorithm.
     221    class BlockSearchPivotRule
     222    {
     223    private:
     224
     225      NetworkSimplex &_ns;
     226      EdgeIt _next_edge, _min_edge;
     227      int _block_size;
     228
     229      static const int MIN_BLOCK_SIZE = 10;
     230
     231    public:
     232
     233      /// Constructor.
     234      BlockSearchPivotRule(NetworkSimplex &ns) :
     235        _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
     236      {
     237        _block_size = 2 * int(sqrt(countEdges(_ns._graph)));
     238        if (_block_size < MIN_BLOCK_SIZE) _block_size = MIN_BLOCK_SIZE;
     239      }
     240
     241      /// Finds the next entering edge.
     242      bool findEnteringEdge() {
     243        Cost curr, min = 0;
     244        int cnt = 0;
     245        for (EdgeIt e = _next_edge; e != INVALID; ++e) {
     246          if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
     247            min = curr;
     248            _min_edge = e;
     249          }
     250          if (++cnt == _block_size) {
     251            if (min < 0) break;
     252            cnt = 0;
     253          }
     254        }
     255        if (min == 0) {
     256          for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
     257            if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
     258              min = curr;
     259              _min_edge = e;
     260            }
     261            if (++cnt == _block_size) {
     262              if (min < 0) break;
     263              cnt = 0;
     264            }
     265          }
     266        }
     267        _ns._in_edge = _min_edge;
     268        _next_edge = ++_min_edge;
     269        return min < 0;
     270      }
     271    }; //class BlockSearchPivotRule
     272
     273    /// \brief Implementation of the "Limited Search" pivot rule for the
     274    /// \ref NetworkSimplex "network simplex" algorithm.
     275    ///
     276    /// This class implements the "Limited Search" pivot rule
     277    /// for the \ref NetworkSimplex "network simplex" algorithm.
     278    class LimitedSearchPivotRule
     279    {
     280    private:
     281
     282      NetworkSimplex &_ns;
     283      EdgeIt _next_edge, _min_edge;
     284      int _sample_size;
     285
     286      static const int MIN_SAMPLE_SIZE = 10;
     287      static const double SAMPLE_SIZE_FACTOR = 0.0015;
     288
     289    public:
     290
     291      /// Constructor.
     292      LimitedSearchPivotRule(NetworkSimplex &ns) :
     293        _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
     294      {
     295        _sample_size = int(SAMPLE_SIZE_FACTOR * countEdges(_ns._graph));
     296        if (_sample_size < MIN_SAMPLE_SIZE)
     297          _sample_size = MIN_SAMPLE_SIZE;
     298      }
     299
     300      /// Finds the next entering edge.
     301      bool findEnteringEdge() {
     302        Cost curr, min = 0;
     303        int cnt = 0;
     304        for (EdgeIt e = _next_edge; e != INVALID; ++e) {
     305          if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
     306            min = curr;
     307            _min_edge = e;
     308          }
     309          if (curr < 0 && ++cnt == _sample_size) break;
     310        }
     311        if (min == 0) {
     312          for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
     313            if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
     314              min = curr;
     315              _min_edge = e;
     316            }
     317            if (curr < 0 && ++cnt == _sample_size) break;
     318          }
     319        }
     320        _ns._in_edge = _min_edge;
     321        _next_edge = ++_min_edge;
     322        return min < 0;
     323      }
     324    }; //class LimitedSearchPivotRule
     325
     326    /// \brief Implementation of the "Candidate List" pivot rule for the
     327    /// \ref NetworkSimplex "network simplex" algorithm.
     328    ///
     329    /// This class implements the "Candidate List" pivot rule
     330    /// for the \ref NetworkSimplex "network simplex" algorithm.
     331    class CandidateListPivotRule
     332    {
     333    private:
     334
     335      NetworkSimplex &_ns;
     336
     337      // The list of candidate edges.
     338      std::vector<Edge> _candidates;
     339      // The maximum length of the edge list.
     340      int _list_length;
     341      // The maximum number of minor iterations between two major
     342      // itarations.
     343      int _minor_limit;
     344
     345      int _minor_count;
     346      EdgeIt _next_edge;
     347
     348      static const double LIST_LENGTH_FACTOR = 0.002;
     349      static const double MINOR_LIMIT_FACTOR = 0.1;
     350      static const int MIN_LIST_LENGTH = 10;
     351      static const int MIN_MINOR_LIMIT = 2;
     352
     353    public:
     354
     355      /// Constructor.
     356      CandidateListPivotRule(NetworkSimplex &ns) :
     357        _ns(ns), _next_edge(ns._graph)
     358      {
     359        int edge_num = countEdges(_ns._graph);
     360        _minor_count = 0;
     361        _list_length = int(edge_num * LIST_LENGTH_FACTOR);
     362        if (_list_length < MIN_LIST_LENGTH)
     363          _list_length = MIN_LIST_LENGTH;
     364        _minor_limit = int(_list_length * MINOR_LIMIT_FACTOR);
     365        if (_minor_limit < MIN_MINOR_LIMIT)
     366          _minor_limit = MIN_MINOR_LIMIT;
     367      }
     368
     369      /// Finds the next entering edge.
     370      bool findEnteringEdge() {
     371        Cost min, curr;
     372        if (_minor_count < _minor_limit && _candidates.size() > 0) {
     373          // Minor iteration
     374          ++_minor_count;
     375          Edge e;
     376          min = 0;
     377          for (int i = 0; i < int(_candidates.size()); ++i) {
     378            e = _candidates[i];
     379            if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
     380              min = curr;
     381              _ns._in_edge = e;
     382            }
     383          }
     384          if (min < 0) return true;
     385        }
     386
     387        // Major iteration
     388        _candidates.clear();
     389        EdgeIt e = _next_edge;
     390        min = 0;
     391        for ( ; e != INVALID; ++e) {
     392          if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
     393            _candidates.push_back(e);
     394            if (curr < min) {
     395              min = curr;
     396              _ns._in_edge = e;
     397            }
     398            if (int(_candidates.size()) == _list_length) break;
     399          }
     400        }
     401        if (int(_candidates.size()) < _list_length) {
     402          for (e = EdgeIt(_ns._graph); e != _next_edge; ++e) {
     403            if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
     404              _candidates.push_back(e);
     405              if (curr < min) {
     406                min = curr;
     407                _ns._in_edge = e;
     408              }
     409              if (int(_candidates.size()) == _list_length) break;
     410            }
     411          }
     412        }
     413        if (_candidates.size() == 0) return false;
     414        _minor_count = 1;
     415        _next_edge = ++e;
     416        return true;
     417      }
     418    }; //class CandidateListPivotRule
     419
     420  private:
     421
     422    // State constant for edges at their lower bounds.
     423    static const int STATE_LOWER =  1;
     424    // State constant for edges in the spanning tree.
     425    static const int STATE_TREE  =  0;
     426    // State constant for edges at their upper bounds.
     427    static const int STATE_UPPER = -1;
     428
     429    // Constant for the combined pivot rule.
     430    static const int COMBINED_PIVOT_MAX_DEG = 5;
     431
     432  private:
     433
     434    // The directed graph the algorithm runs on
     435    SGraph _graph;
     436    // The original graph
     437    const Graph &_graph_ref;
     438    // The original lower bound map
     439    const LowerMap *_lower;
     440    // The capacity map
     441    SCapacityMap _capacity;
     442    // The cost map
     443    SCostMap _cost;
     444    // The supply map
     445    SSupplyMap _supply;
     446    bool _valid_supply;
     447
     448    // Edge map of the current flow
     449    SCapacityMap _flow;
     450    // Node map of the current potentials
     451    SPotentialMap _potential;
     452
     453    // The depth node map of the spanning tree structure
     454    IntNodeMap _depth;
     455    // The parent node map of the spanning tree structure
     456    NodeNodeMap _parent;
     457    // The pred_edge node map of the spanning tree structure
     458    EdgeNodeMap _pred_edge;
     459    // The thread node map of the spanning tree structure
     460    NodeNodeMap _thread;
     461    // The forward node map of the spanning tree structure
     462    BoolNodeMap _forward;
     463    // The state edge map
     464    IntEdgeMap _state;
     465    // The root node of the starting spanning tree
     466    Node _root;
     467
     468    // The reduced cost map
     469    ReducedCostMap _red_cost;
     470
     471    // Members for handling the original graph
     472    FlowMap _flow_result;
     473    PotentialMap _potential_result;
     474    NodeRefMap _node_ref;
     475    EdgeRefMap _edge_ref;
    200476
    201477    // The entering edge of the current pivot iteration.
    202     Edge in_edge;
     478    Edge _in_edge;
     479
    203480    // Temporary nodes used in the current pivot iteration.
    204481    Node join, u_in, v_in, u_out, v_out;
     
    209486    Capacity delta;
    210487
    211 #ifdef EDGE_BLOCK_PIVOT
    212     /// The size of blocks for the "Edge Block" pivot rule.
    213     int block_size;
    214 #endif
    215 #ifdef CANDIDATE_LIST_PIVOT
    216     /// \brief The list of candidate edges for the "Candidate List"
    217     /// pivot rule.
    218     std::vector<Edge> candidates;
    219     /// \brief The maximum length of the edge list for the
    220     /// "Candidate List" pivot rule.
    221     int list_length;
    222     /// \brief The maximum number of minor iterations between two major
    223     /// itarations.
    224     int minor_limit;
    225     /// \brief The number of minor iterations.
    226     int minor_count;
    227 #endif
    228 #ifdef SORTED_LIST_PIVOT
    229     /// \brief The list of candidate edges for the "Sorted List"
    230     /// pivot rule.
    231     std::vector<Edge> candidates;
    232     /// \brief The maximum length of the edge list for the
    233     /// "Sorted List" pivot rule.
    234     int list_length;
    235     int list_index;
    236 #endif
    237 
    238488  public :
    239489
     
    242492    /// General constructor of the class (with lower bounds).
    243493    ///
    244     /// \param _graph The directed graph the algorithm runs on.
    245     /// \param _lower The lower bounds of the edges.
    246     /// \param _capacity The capacities (upper bounds) of the edges.
    247     /// \param _cost The cost (length) values of the edges.
    248     /// \param _supply The supply values of the nodes (signed).
    249     NetworkSimplex( const Graph &_graph,
    250                     const LowerMap &_lower,
    251                     const CapacityMap &_capacity,
    252                     const CostMap &_cost,
    253                     const SupplyMap &_supply ) :
    254       graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),
    255       supply(graph), flow(graph), flow_result(_graph), potential(graph),
    256       potential_result(_graph), depth(graph), parent(graph),
    257       pred_edge(graph), thread(graph), forward(graph), state(graph),
    258       node_ref(graph_ref), edge_ref(graph_ref),
    259       red_cost(graph, cost, potential)
     494    /// \param graph The directed graph the algorithm runs on.
     495    /// \param lower The lower bounds of the edges.
     496    /// \param capacity The capacities (upper bounds) of the edges.
     497    /// \param cost The cost (length) values of the edges.
     498    /// \param supply The supply values of the nodes (signed).
     499    NetworkSimplex( const Graph &graph,
     500                    const LowerMap &lower,
     501                    const CapacityMap &capacity,
     502                    const CostMap &cost,
     503                    const SupplyMap &supply ) :
     504      _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
     505      _cost(_graph), _supply(_graph), _flow(_graph),
     506      _potential(_graph), _depth(_graph), _parent(_graph),
     507      _pred_edge(_graph), _thread(_graph), _forward(_graph),
     508      _state(_graph), _red_cost(_graph, _cost, _potential),
     509      _flow_result(graph), _potential_result(graph),
     510      _node_ref(graph), _edge_ref(graph)
    260511    {
    261512      // Checking the sum of supply values
    262513      Supply sum = 0;
    263       for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
    264         sum += _supply[n];
    265       if (!(valid_supply = sum == 0)) return;
    266 
    267       // Copying graph_ref to graph
    268       graph.reserveNode(countNodes(graph_ref) + 1);
    269       graph.reserveEdge(countEdges(graph_ref) + countNodes(graph_ref));
    270       copyGraph(graph, graph_ref)
    271         .edgeMap(cost, _cost)
    272         .nodeRef(node_ref)
    273         .edgeRef(edge_ref)
     514      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
     515        sum += supply[n];
     516      if (!(_valid_supply = sum == 0)) return;
     517
     518      // Copying _graph_ref to _graph
     519      _graph.reserveNode(countNodes(_graph_ref) + 1);
     520      _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
     521      copyGraph(_graph, _graph_ref)
     522        .edgeMap(_cost, cost)
     523        .nodeRef(_node_ref)
     524        .edgeRef(_edge_ref)
    274525        .run();
    275526
    276527      // Removing non-zero lower bounds
    277       for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {
    278         capacity[edge_ref[e]] = _capacity[e] - _lower[e];
    279       }
    280       for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {
    281         Supply s = _supply[n];
    282         for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
    283           s += _lower[e];
    284         for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
    285           s -= _lower[e];
    286         supply[node_ref[n]] = s;
     528      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
     529        _capacity[_edge_ref[e]] = capacity[e] - lower[e];
     530      }
     531      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {
     532        Supply s = supply[n];
     533        for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)
     534          s += lower[e];
     535        for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)
     536          s -= lower[e];
     537        _supply[_node_ref[n]] = s;
    287538      }
    288539    }
     
    292543    /// General constructor of the class (without lower bounds).
    293544    ///
    294     /// \param _graph The directed graph the algorithm runs on.
    295     /// \param _capacity The capacities (upper bounds) of the edges.
    296     /// \param _cost The cost (length) values of the edges.
    297     /// \param _supply The supply values of the nodes (signed).
    298     NetworkSimplex( const Graph &_graph,
    299                     const CapacityMap &_capacity,
    300                     const CostMap &_cost,
    301                     const SupplyMap &_supply ) :
    302       graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),
    303       supply(graph), flow(graph), flow_result(_graph), potential(graph),
    304       potential_result(_graph), depth(graph), parent(graph),
    305       pred_edge(graph), thread(graph), forward(graph), state(graph),
    306       node_ref(graph_ref), edge_ref(graph_ref),
    307       red_cost(graph, cost, potential)
     545    /// \param graph The directed graph the algorithm runs on.
     546    /// \param capacity The capacities (upper bounds) of the edges.
     547    /// \param cost The cost (length) values of the edges.
     548    /// \param supply The supply values of the nodes (signed).
     549    NetworkSimplex( const Graph &graph,
     550                    const CapacityMap &capacity,
     551                    const CostMap &cost,
     552                    const SupplyMap &supply ) :
     553      _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
     554      _cost(_graph), _supply(_graph), _flow(_graph),
     555      _potential(_graph), _depth(_graph), _parent(_graph),
     556      _pred_edge(_graph), _thread(_graph), _forward(_graph),
     557      _state(_graph), _red_cost(_graph, _cost, _potential),
     558      _flow_result(graph), _potential_result(graph),
     559      _node_ref(graph), _edge_ref(graph)
    308560    {
    309561      // Checking the sum of supply values
    310562      Supply sum = 0;
    311       for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
    312         sum += _supply[n];
    313       if (!(valid_supply = sum == 0)) return;
    314 
    315       // Copying graph_ref to graph
    316       copyGraph(graph, graph_ref)
    317         .edgeMap(capacity, _capacity)
    318         .edgeMap(cost, _cost)
    319         .nodeMap(supply, _supply)
    320         .nodeRef(node_ref)
    321         .edgeRef(edge_ref)
     563      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
     564        sum += supply[n];
     565      if (!(_valid_supply = sum == 0)) return;
     566
     567      // Copying _graph_ref to graph
     568      copyGraph(_graph, _graph_ref)
     569        .edgeMap(_capacity, capacity)
     570        .edgeMap(_cost, cost)
     571        .nodeMap(_supply, supply)
     572        .nodeRef(_node_ref)
     573        .edgeRef(_edge_ref)
    322574        .run();
    323575    }
     
    327579    /// Simple constructor of the class (with lower bounds).
    328580    ///
    329     /// \param _graph The directed graph the algorithm runs on.
    330     /// \param _lower The lower bounds of the edges.
    331     /// \param _capacity The capacities (upper bounds) of the edges.
    332     /// \param _cost The cost (length) values of the edges.
    333     /// \param _s The source node.
    334     /// \param _t The target node.
    335     /// \param _flow_value The required amount of flow from node \c _s
    336     /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t).
    337     NetworkSimplex( const Graph &_graph,
    338                     const LowerMap &_lower,
    339                     const CapacityMap &_capacity,
    340                     const CostMap &_cost,
    341                     typename Graph::Node _s,
    342                     typename Graph::Node _t,
    343                     typename SupplyMap::Value _flow_value ) :
    344       graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),
    345       supply(graph), flow(graph), flow_result(_graph), potential(graph),
    346       potential_result(_graph), depth(graph), parent(graph),
    347       pred_edge(graph), thread(graph), forward(graph), state(graph),
    348       node_ref(graph_ref), edge_ref(graph_ref),
    349       red_cost(graph, cost, potential)
     581    /// \param graph The directed graph the algorithm runs on.
     582    /// \param lower The lower bounds of the edges.
     583    /// \param capacity The capacities (upper bounds) of the edges.
     584    /// \param cost The cost (length) values of the edges.
     585    /// \param s The source node.
     586    /// \param t The target node.
     587    /// \param flow_value The required amount of flow from node \c s
     588    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
     589    NetworkSimplex( const Graph &graph,
     590                    const LowerMap &lower,
     591                    const CapacityMap &capacity,
     592                    const CostMap &cost,
     593                    typename Graph::Node s,
     594                    typename Graph::Node t,
     595                    typename SupplyMap::Value flow_value ) :
     596      _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
     597      _cost(_graph), _supply(_graph), _flow(_graph),
     598      _potential(_graph), _depth(_graph), _parent(_graph),
     599      _pred_edge(_graph), _thread(_graph), _forward(_graph),
     600      _state(_graph), _red_cost(_graph, _cost, _potential),
     601      _flow_result(graph), _potential_result(graph),
     602      _node_ref(graph), _edge_ref(graph)
    350603    {
    351       // Copying graph_ref to graph
    352       copyGraph(graph, graph_ref)
    353         .edgeMap(cost, _cost)
    354         .nodeRef(node_ref)
    355         .edgeRef(edge_ref)
     604      // Copying _graph_ref to graph
     605      copyGraph(_graph, _graph_ref)
     606        .edgeMap(_cost, cost)
     607        .nodeRef(_node_ref)
     608        .edgeRef(_edge_ref)
    356609        .run();
    357610
    358611      // Removing non-zero lower bounds
    359       for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {
    360         capacity[edge_ref[e]] = _capacity[e] - _lower[e];
    361       }
    362       for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {
    363         Supply s = 0;
    364         if (n == _s) s =  _flow_value;
    365         if (n == _t) s = -_flow_value;
    366         for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
    367           s += _lower[e];
    368         for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
    369           s -= _lower[e];
    370         supply[node_ref[n]] = s;
    371       }
    372       valid_supply = true;
     612      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
     613        _capacity[_edge_ref[e]] = capacity[e] - lower[e];
     614      }
     615      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {
     616        Supply sum = 0;
     617        if (n == s) sum =  flow_value;
     618        if (n == t) sum = -flow_value;
     619        for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)
     620          sum += lower[e];
     621        for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)
     622          sum -= lower[e];
     623        _supply[_node_ref[n]] = sum;
     624      }
     625      _valid_supply = true;
    373626    }
    374627
     
    377630    /// Simple constructor of the class (without lower bounds).
    378631    ///
    379     /// \param _graph The directed graph the algorithm runs on.
    380     /// \param _capacity The capacities (upper bounds) of the edges.
    381     /// \param _cost The cost (length) values of the edges.
    382     /// \param _s The source node.
    383     /// \param _t The target node.
    384     /// \param _flow_value The required amount of flow from node \c _s
    385     /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t).
    386     NetworkSimplex( const Graph &_graph,
    387                     const CapacityMap &_capacity,
    388                     const CostMap &_cost,
    389                     typename Graph::Node _s,
    390                     typename Graph::Node _t,
    391                     typename SupplyMap::Value _flow_value ) :
    392       graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),
    393       supply(graph, 0), flow(graph), flow_result(_graph), potential(graph),
    394       potential_result(_graph), depth(graph), parent(graph),
    395       pred_edge(graph), thread(graph), forward(graph), state(graph),
    396       node_ref(graph_ref), edge_ref(graph_ref),
    397       red_cost(graph, cost, potential)
     632    /// \param graph The directed graph the algorithm runs on.
     633    /// \param capacity The capacities (upper bounds) of the edges.
     634    /// \param cost The cost (length) values of the edges.
     635    /// \param s The source node.
     636    /// \param t The target node.
     637    /// \param flow_value The required amount of flow from node \c s
     638    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
     639    NetworkSimplex( const Graph &graph,
     640                    const CapacityMap &capacity,
     641                    const CostMap &cost,
     642                    typename Graph::Node s,
     643                    typename Graph::Node t,
     644                    typename SupplyMap::Value flow_value ) :
     645      _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
     646      _cost(_graph), _supply(_graph, 0), _flow(_graph),
     647      _potential(_graph), _depth(_graph), _parent(_graph),
     648      _pred_edge(_graph), _thread(_graph), _forward(_graph),
     649      _state(_graph), _red_cost(_graph, _cost, _potential),
     650      _flow_result(graph), _potential_result(graph),
     651      _node_ref(graph), _edge_ref(graph)
    398652    {
    399       // Copying graph_ref to graph
    400       copyGraph(graph, graph_ref)
    401         .edgeMap(capacity, _capacity)
    402         .edgeMap(cost, _cost)
    403         .nodeRef(node_ref)
    404         .edgeRef(edge_ref)
     653      // Copying _graph_ref to graph
     654      copyGraph(_graph, _graph_ref)
     655        .edgeMap(_capacity, capacity)
     656        .edgeMap(_cost, cost)
     657        .nodeRef(_node_ref)
     658        .edgeRef(_edge_ref)
    405659        .run();
    406       supply[node_ref[_s]] =  _flow_value;
    407       supply[node_ref[_t]] = -_flow_value;
    408       valid_supply = true;
     660      _supply[_node_ref[s]] =  flow_value;
     661      _supply[_node_ref[t]] = -flow_value;
     662      _valid_supply = true;
    409663    }
    410664
     
    413667    /// Runs the algorithm.
    414668    ///
     669    /// \param pivot_rule The pivot rule that is used during the
     670    /// algorithm.
     671    ///
     672    /// The available pivot rules:
     673    ///
     674    /// - FIRST_ELIGIBLE_PIVOT The next eligible edge is selected in
     675    /// a wraparound fashion in every iteration
     676    /// (\ref FirstEligiblePivotRule).
     677    ///
     678    /// - BEST_ELIGIBLE_PIVOT The best eligible edge is selected in
     679    /// every iteration (\ref BestEligiblePivotRule).
     680    ///
     681    /// - BLOCK_SEARCH_PIVOT A specified number of edges are examined in
     682    /// every iteration in a wraparound fashion and the best eligible
     683    /// edge is selected from this block (\ref BlockSearchPivotRule).
     684    ///
     685    /// - LIMITED_SEARCH_PIVOT A specified number of eligible edges are
     686    /// examined in every iteration in a wraparound fashion and the best
     687    /// one is selected from them (\ref LimitedSearchPivotRule).
     688    ///
     689    /// - CANDIDATE_LIST_PIVOT In major iterations a candidate list is
     690    /// built from eligible edges and it is used for edge selection in
     691    /// the following minor iterations (\ref CandidateListPivotRule).
     692    ///
     693    /// - COMBINED_PIVOT This is a combined version of the two fastest
     694    /// pivot rules.
     695    /// For rather sparse graphs \ref LimitedSearchPivotRule
     696    /// "Limited Search" implementation is used, otherwise
     697    /// \ref BlockSearchPivotRule "Block Search" pivot rule is used.
     698    /// According to our benchmark tests this combined method is the
     699    /// most efficient.
     700    ///
    415701    /// \return \c true if a feasible flow can be found.
    416     bool run() {
    417       return init() && start();
    418     }
    419 
    420     /// \brief Returns a const reference to the flow map.
    421     ///
    422     /// Returns a const reference to the flow map.
     702    bool run(PivotRuleEnum pivot_rule = COMBINED_PIVOT) {
     703      return init() && start(pivot_rule);
     704    }
     705
     706    /// \brief Returns a const reference to the edge map storing the
     707    /// found flow.
     708    ///
     709    /// Returns a const reference to the edge map storing the found flow.
    423710    ///
    424711    /// \pre \ref run() must be called before using this function.
    425712    const FlowMap& flowMap() const {
    426       return flow_result;
    427     }
    428 
    429     /// \brief Returns a const reference to the potential map (the dual
    430     /// solution).
    431     ///
    432     /// Returns a const reference to the potential map (the dual
    433     /// solution).
     713      return _flow_result;
     714    }
     715
     716    /// \brief Returns a const reference to the node map storing the
     717    /// found potentials (the dual solution).
     718    ///
     719    /// Returns a const reference to the node map storing the found
     720    /// potentials (the dual solution).
    434721    ///
    435722    /// \pre \ref run() must be called before using this function.
    436723    const PotentialMap& potentialMap() const {
    437       return potential_result;
     724      return _potential_result;
    438725    }
    439726
     
    446733    Cost totalCost() const {
    447734      Cost c = 0;
    448       for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
    449         c += flow_result[e] * cost[edge_ref[e]];
     735      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
     736        c += _flow_result[e] * _cost[_edge_ref[e]];
    450737      return c;
    451738    }
    452739
    453   protected:
     740  private:
    454741
    455742    /// \brief Extends the underlaying graph and initializes all the
    456743    /// node and edge maps.
    457744    bool init() {
    458       if (!valid_supply) return false;
     745      if (!_valid_supply) return false;
    459746
    460747      // Initializing state and flow maps
    461       for (EdgeIt e(graph); e != INVALID; ++e) {
    462         flow[e] = 0;
    463         state[e] = LOWER;
     748      for (EdgeIt e(_graph); e != INVALID; ++e) {
     749        _flow[e] = 0;
     750        _state[e] = STATE_LOWER;
    464751      }
    465752
    466753      // Adding an artificial root node to the graph
    467       root = graph.addNode();
    468       parent[root] = INVALID;
    469       pred_edge[root] = INVALID;
    470       depth[root] = 0;
    471       supply[root] = 0;
    472       potential[root] = 0;
     754      _root = _graph.addNode();
     755      _parent[_root] = INVALID;
     756      _pred_edge[_root] = INVALID;
     757      _depth[_root] = 0;
     758      _supply[_root] = 0;
     759      _potential[_root] = 0;
    473760
    474761      // Adding artificial edges to the graph and initializing the node
    475762      // maps of the spanning tree data structure
    476       Supply sum = 0;
    477       Node last = root;
     763      Node last = _root;
    478764      Edge e;
    479765      Cost max_cost = std::numeric_limits<Cost>::max() / 4;
    480       for (NodeIt u(graph); u != INVALID; ++u) {
    481         if (u == root) continue;
    482         thread[last] = u;
     766      for (NodeIt u(_graph); u != INVALID; ++u) {
     767        if (u == _root) continue;
     768        _thread[last] = u;
    483769        last = u;
    484         parent[u] = root;
    485         depth[u] = 1;
    486         sum += supply[u];
    487         if (supply[u] >= 0) {
    488           e = graph.addEdge(u, root);
    489           flow[e] = supply[u];
    490           forward[u] = true;
    491           potential[u] = max_cost;
     770        _parent[u] = _root;
     771        _depth[u] = 1;
     772        if (_supply[u] >= 0) {
     773          e = _graph.addEdge(u, _root);
     774          _flow[e] = _supply[u];
     775          _forward[u] = true;
     776          _potential[u] = -max_cost;
    492777        } else {
    493           e = graph.addEdge(root, u);
    494           flow[e] = -supply[u];
    495           forward[u] = false;
    496           potential[u] = -max_cost;
    497         }
    498         cost[e] = max_cost;
    499         capacity[e] = std::numeric_limits<Capacity>::max();
    500         state[e] = TREE;
    501         pred_edge[u] = e;
    502       }
    503       thread[last] = root;
    504 
    505 #ifdef EDGE_BLOCK_PIVOT
    506       // Initializing block_size for the edge block pivot rule
    507       int edge_num = countEdges(graph);
    508       block_size = 2 * int(sqrt(countEdges(graph)));
    509       if (block_size < MIN_BLOCK_SIZE) block_size = MIN_BLOCK_SIZE;
    510 #endif
    511 #ifdef CANDIDATE_LIST_PIVOT
    512       int edge_num = countEdges(graph);
    513       minor_count = 0;
    514       list_length = edge_num / LIST_LENGTH_DIV;
    515       minor_limit = edge_num / MINOR_LIMIT_DIV;
    516 #endif
    517 #ifdef SORTED_LIST_PIVOT
    518       int edge_num = countEdges(graph);
    519       list_index = 0;
    520       list_length = edge_num / LIST_LENGTH_DIV;
    521 #endif
    522 
    523       return sum == 0;
    524     }
    525 
    526 #ifdef FIRST_ELIGIBLE_PIVOT
    527     /// \brief Finds entering edge according to the "First Eligible"
    528     /// pivot rule.
    529     bool findEnteringEdge(EdgeIt &next_edge) {
    530       for (EdgeIt e = next_edge; e != INVALID; ++e) {
    531         if (state[e] * red_cost[e] < 0) {
    532           in_edge = e;
    533           next_edge = ++e;
    534           return true;
    535         }
    536       }
    537       for (EdgeIt e(graph); e != next_edge; ++e) {
    538         if (state[e] * red_cost[e] < 0) {
    539           in_edge = e;
    540           next_edge = ++e;
    541           return true;
    542         }
    543       }
    544       return false;
    545     }
    546 #endif
    547 
    548 #ifdef BEST_ELIGIBLE_PIVOT
    549     /// \brief Finds entering edge according to the "Best Eligible"
    550     /// pivot rule.
    551     bool findEnteringEdge() {
    552       Cost min = 0;
    553       for (EdgeIt e(graph); e != INVALID; ++e) {
    554         if (state[e] * red_cost[e] < min) {
    555           min = state[e] * red_cost[e];
    556           in_edge = e;
    557         }
    558       }
    559       return min < 0;
    560     }
    561 #endif
    562 
    563 #ifdef EDGE_BLOCK_PIVOT
    564     /// \brief Finds entering edge according to the "Edge Block"
    565     /// pivot rule.
    566     bool findEnteringEdge(EdgeIt &next_edge) {
    567       // Performing edge block selection
    568       Cost curr, min = 0;
    569       EdgeIt min_edge(graph);
    570       int cnt = 0;
    571       for (EdgeIt e = next_edge; e != INVALID; ++e) {
    572         if ((curr = state[e] * red_cost[e]) < min) {
    573           min = curr;
    574           min_edge = e;
    575         }
    576         if (++cnt == block_size) {
    577           if (min < 0) break;
    578           cnt = 0;
    579         }
    580       }
    581       if (!(min < 0)) {
    582         for (EdgeIt e(graph); e != next_edge; ++e) {
    583           if ((curr = state[e] * red_cost[e]) < min) {
    584             min = curr;
    585             min_edge = e;
    586           }
    587           if (++cnt == block_size) {
    588             if (min < 0) break;
    589             cnt = 0;
    590           }
    591         }
    592       }
    593       in_edge = min_edge;
    594       if ((next_edge = ++min_edge) == INVALID)
    595         next_edge = EdgeIt(graph);
    596       return min < 0;
    597     }
    598 #endif
    599 
    600 #ifdef CANDIDATE_LIST_PIVOT
    601     /// \brief Finds entering edge according to the "Candidate List"
    602     /// pivot rule.
    603     bool findEnteringEdge() {
    604       typedef typename std::vector<Edge>::iterator ListIt;
    605 
    606       if (minor_count >= minor_limit || candidates.size() == 0) {
    607         // Major iteration
    608         candidates.clear();
    609         for (EdgeIt e(graph); e != INVALID; ++e) {
    610           if (state[e] * red_cost[e] < 0) {
    611             candidates.push_back(e);
    612             if (candidates.size() == list_length) break;
    613           }
    614         }
    615         if (candidates.size() == 0) return false;
    616       }
    617 
    618       // Minor iteration
    619       ++minor_count;
    620       Cost min = 0;
    621       Edge e;
    622       for (int i = 0; i < candidates.size(); ++i) {
    623         e = candidates[i];
    624         if (state[e] * red_cost[e] < min) {
    625           min = state[e] * red_cost[e];
    626           in_edge = e;
    627         }
    628       }
     778          e = _graph.addEdge(_root, u);
     779          _flow[e] = -_supply[u];
     780          _forward[u] = false;
     781          _potential[u] = max_cost;
     782        }
     783        _cost[e] = max_cost;
     784        _capacity[e] = std::numeric_limits<Capacity>::max();
     785        _state[e] = STATE_TREE;
     786        _pred_edge[u] = e;
     787      }
     788      _thread[last] = _root;
     789
    629790      return true;
    630791    }
    631 #endif
    632 
    633 #ifdef SORTED_LIST_PIVOT
    634     /// \brief Functor class to compare edges during sort of the
    635     /// candidate list.
    636     class SortFunc
    637     {
    638     private:
    639       const IntEdgeMap &st;
    640       const ReducedCostMap &rc;
    641     public:
    642       SortFunc(const IntEdgeMap &_st, const ReducedCostMap &_rc) :
    643         st(_st), rc(_rc) {}
    644       bool operator()(const Edge &e1, const Edge &e2) {
    645         return st[e1] * rc[e1] < st[e2] * rc[e2];
    646       }
    647     };
    648 
    649     /// \brief Finds entering edge according to the "Sorted List"
    650     /// pivot rule.
    651     bool findEnteringEdge() {
    652       static SortFunc sort_func(state, red_cost);
    653 
    654       // Minor iteration
    655       while (list_index < candidates.size()) {
    656         in_edge = candidates[list_index++];
    657         if (state[in_edge] * red_cost[in_edge] < 0) return true;
    658       }
    659 
    660       // Major iteration
    661       candidates.clear();
    662       Cost curr, min = 0;
    663       for (EdgeIt e(graph); e != INVALID; ++e) {
    664         if ((curr = state[e] * red_cost[e]) < min / LOWER_DIV) {
    665           candidates.push_back(e);
    666           if (curr < min) min = curr;
    667           if (candidates.size() == list_length) break;
    668         }
    669       }
    670       if (candidates.size() == 0) return false;
    671       sort(candidates.begin(), candidates.end(), sort_func);
    672       in_edge = candidates[0];
    673       list_index = 1;
    674       return true;
    675     }
    676 #endif
    677 
    678     /// \brief Finds the join node.
     792
     793    /// Finds the join node.
    679794    Node findJoinNode() {
    680       // Finding the join node
    681       Node u = graph.source(in_edge);
    682       Node v = graph.target(in_edge);
     795      Node u = _graph.source(_in_edge);
     796      Node v = _graph.target(_in_edge);
    683797      while (u != v) {
    684         if (depth[u] == depth[v]) {
    685           u = parent[u];
    686           v = parent[v];
    687         }
    688         else if (depth[u] > depth[v]) u = parent[u];
    689         else v = parent[v];
     798        if (_depth[u] == _depth[v]) {
     799          u = _parent[u];
     800          v = _parent[v];
     801        }
     802        else if (_depth[u] > _depth[v]) u = _parent[u];
     803        else v = _parent[v];
    690804      }
    691805      return u;
     
    697811      // Initializing first and second nodes according to the direction
    698812      // of the cycle
    699       if (state[in_edge] == LOWER) {
    700         first = graph.source(in_edge);
    701         second  = graph.target(in_edge);
     813      if (_state[_in_edge] == STATE_LOWER) {
     814        first  = _graph.source(_in_edge);
     815        second = _graph.target(_in_edge);
    702816      } else {
    703         first = graph.target(in_edge);
    704         second  = graph.source(in_edge);
    705       }
    706       delta = capacity[in_edge];
     817        first  = _graph.target(_in_edge);
     818        second = _graph.source(_in_edge);
     819      }
     820      delta = _capacity[_in_edge];
    707821      bool result = false;
    708822      Capacity d;
     
    711825      // Searching the cycle along the path form the first node to the
    712826      // root node
    713       for (Node u = first; u != join; u = parent[u]) {
    714         e = pred_edge[u];
    715         d = forward[u] ? flow[e] : capacity[e] - flow[e];
     827      for (Node u = first; u != join; u = _parent[u]) {
     828        e = _pred_edge[u];
     829        d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e];
    716830        if (d < delta) {
    717831          delta = d;
     
    724838      // Searching the cycle along the path form the second node to the
    725839      // root node
    726       for (Node u = second; u != join; u = parent[u]) {
    727         e = pred_edge[u];
    728         d = forward[u] ? capacity[e] - flow[e] : flow[e];
     840      for (Node u = second; u != join; u = _parent[u]) {
     841        e = _pred_edge[u];
     842        d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e];
    729843        if (d <= delta) {
    730844          delta = d;
     
    738852    }
    739853
    740     /// \brief Changes \c flow and \c state edge maps.
     854    /// Changes \c flow and \c state edge maps.
    741855    void changeFlows(bool change) {
    742856      // Augmenting along the cycle
    743857      if (delta > 0) {
    744         Capacity val = state[in_edge] * delta;
    745         flow[in_edge] += val;
    746         for (Node u = graph.source(in_edge); u != join; u = parent[u]) {
    747           flow[pred_edge[u]] += forward[u] ? -val : val;
    748         }
    749         for (Node u = graph.target(in_edge); u != join; u = parent[u]) {
    750           flow[pred_edge[u]] += forward[u] ? val : -val;
     858        Capacity val = _state[_in_edge] * delta;
     859        _flow[_in_edge] += val;
     860        for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) {
     861          _flow[_pred_edge[u]] += _forward[u] ? -val : val;
     862        }
     863        for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) {
     864          _flow[_pred_edge[u]] += _forward[u] ? val : -val;
    751865        }
    752866      }
    753867      // Updating the state of the entering and leaving edges
    754868      if (change) {
    755         state[in_edge] = TREE;
    756         state[pred_edge[u_out]] =
    757           (flow[pred_edge[u_out]] == 0) ? LOWER : UPPER;
     869        _state[_in_edge] = STATE_TREE;
     870        _state[_pred_edge[u_out]] =
     871          (_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
    758872      } else {
    759         state[in_edge] = -state[in_edge];
    760       }
    761     }
    762 
    763     /// \brief Updates \c thread and \c parent node maps.
     873        _state[_in_edge] = -_state[_in_edge];
     874      }
     875    }
     876
     877    /// Updates \c thread and \c parent node maps.
    764878    void updateThreadParent() {
    765879      Node u;
    766       v_out = parent[u_out];
     880      v_out = _parent[u_out];
    767881
    768882      // Handling the case when join and v_out coincide
    769883      bool par_first = false;
    770884      if (join == v_out) {
    771         for (u = join; u != u_in && u != v_in; u = thread[u]) ;
     885        for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
    772886        if (u == v_in) {
    773887          par_first = true;
    774           while (thread[u] != u_out) u = thread[u];
     888          while (_thread[u] != u_out) u = _thread[u];
    775889          first = u;
    776890        }
     
    779893      // Finding the last successor of u_in (u) and the node after it
    780894      // (right) according to the thread index
    781       for (u = u_in; depth[thread[u]] > depth[u_in]; u = thread[u]) ;
    782       right = thread[u];
    783       if (thread[v_in] == u_out) {
    784         for (last = u; depth[last] > depth[u_out]; last = thread[last]) ;
    785         if (last == u_out) last = thread[last];
    786       }
    787       else last = thread[v_in];
     895      for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
     896      right = _thread[u];
     897      if (_thread[v_in] == u_out) {
     898        for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;
     899        if (last == u_out) last = _thread[last];
     900      }
     901      else last = _thread[v_in];
    788902
    789903      // Updating stem nodes
    790       thread[v_in] = stem = u_in;
     904      _thread[v_in] = stem = u_in;
    791905      par_stem = v_in;
    792906      while (stem != u_out) {
    793         thread[u] = new_stem = parent[stem];
     907        _thread[u] = new_stem = _parent[stem];
    794908
    795909        // Finding the node just before the stem node (u) according to
    796910        // the original thread index
    797         for (u = new_stem; thread[u] != stem; u = thread[u]) ;
    798         thread[u] = right;
     911        for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
     912        _thread[u] = right;
    799913
    800914        // Changing the parent node of stem and shifting stem and
    801915        // par_stem nodes
    802         parent[stem] = par_stem;
     916        _parent[stem] = par_stem;
    803917        par_stem = stem;
    804918        stem = new_stem;
     
    806920        // Finding the last successor of stem (u) and the node after it
    807921        // (right) according to the thread index
    808         for (u = stem; depth[thread[u]] > depth[stem]; u = thread[u]) ;
    809         right = thread[u];
    810       }
    811       parent[u_out] = par_stem;
    812       thread[u] = last;
     922        for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
     923        right = _thread[u];
     924      }
     925      _parent[u_out] = par_stem;
     926      _thread[u] = last;
    813927
    814928      if (join == v_out && par_first) {
    815         if (first != v_in) thread[first] = right;
     929        if (first != v_in) _thread[first] = right;
    816930      } else {
    817         for (u = v_out; thread[u] != u_out; u = thread[u]) ;
    818         thread[u] = right;
    819       }
    820     }
    821 
    822     /// \brief Updates \c pred_edge and \c forward node maps.
     931        for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
     932        _thread[u] = right;
     933      }
     934    }
     935
     936    /// Updates \c pred_edge and \c forward node maps.
    823937    void updatePredEdge() {
    824938      Node u = u_out, v;
    825939      while (u != u_in) {
    826         v = parent[u];
    827         pred_edge[u] = pred_edge[v];
    828         forward[u] = !forward[v];
     940        v = _parent[u];
     941        _pred_edge[u] = _pred_edge[v];
     942        _forward[u] = !_forward[v];
    829943        u = v;
    830944      }
    831       pred_edge[u_in] = in_edge;
    832       forward[u_in] = (u_in == graph.source(in_edge));
    833     }
    834 
    835     /// \brief Updates \c depth and \c potential node maps.
     945      _pred_edge[u_in] = _in_edge;
     946      _forward[u_in] = (u_in == _graph.source(_in_edge));
     947    }
     948
     949    /// Updates \c depth and \c potential node maps.
    836950    void updateDepthPotential() {
    837       depth[u_in] = depth[v_in] + 1;
    838       potential[u_in] = forward[u_in] ?
    839         potential[v_in] + cost[pred_edge[u_in]] :
    840         potential[v_in] - cost[pred_edge[u_in]];
    841 
    842       Node u = thread[u_in], v;
     951      _depth[u_in] = _depth[v_in] + 1;
     952      _potential[u_in] = _forward[u_in] ?
     953        _potential[v_in] - _cost[_pred_edge[u_in]] :
     954        _potential[v_in] + _cost[_pred_edge[u_in]];
     955
     956      Node u = _thread[u_in], v;
    843957      while (true) {
    844         v = parent[u];
     958        v = _parent[u];
    845959        if (v == INVALID) break;
    846         depth[u] = depth[v] + 1;
    847         potential[u] = forward[u] ?
    848           potential[v] + cost[pred_edge[u]] :
    849           potential[v] - cost[pred_edge[u]];
    850         if (depth[u] <= depth[v_in]) break;
    851         u = thread[u];
    852       }
    853     }
    854 
    855     /// \brief Executes the algorithm.
     960        _depth[u] = _depth[v] + 1;
     961        _potential[u] = _forward[u] ?
     962          _potential[v] - _cost[_pred_edge[u]] :
     963          _potential[v] + _cost[_pred_edge[u]];
     964        if (_depth[u] <= _depth[v_in]) break;
     965        u = _thread[u];
     966      }
     967    }
     968
     969    /// Executes the algorithm.
     970    bool start(PivotRuleEnum pivot_rule) {
     971      switch (pivot_rule) {
     972        case FIRST_ELIGIBLE_PIVOT:
     973          return start<FirstEligiblePivotRule>();
     974        case BEST_ELIGIBLE_PIVOT:
     975          return start<BestEligiblePivotRule>();
     976        case BLOCK_SEARCH_PIVOT:
     977          return start<BlockSearchPivotRule>();
     978        case LIMITED_SEARCH_PIVOT:
     979          return start<LimitedSearchPivotRule>();
     980        case CANDIDATE_LIST_PIVOT:
     981          return start<CandidateListPivotRule>();
     982        case COMBINED_PIVOT:
     983          if ( countEdges(_graph) / countNodes(_graph) <=
     984               COMBINED_PIVOT_MAX_DEG )
     985            return start<LimitedSearchPivotRule>();
     986          else
     987            return start<BlockSearchPivotRule>();
     988      }
     989      return false;
     990    }
     991
     992    template<class PivotRuleImplementation>
    856993    bool start() {
    857       // Processing pivots
    858 #ifdef _DEBUG_ITER_
    859       int iter_num = 0;
    860 #endif
    861 #if defined(FIRST_ELIGIBLE_PIVOT) || defined(EDGE_BLOCK_PIVOT)
    862       EdgeIt next_edge(graph);
    863       while (findEnteringEdge(next_edge))
    864 #else
    865       while (findEnteringEdge())
    866 #endif
    867       {
     994      PivotRuleImplementation pivot(*this);
     995
     996      // Executing the network simplex algorithm
     997      while (pivot.findEnteringEdge()) {
    868998        join = findJoinNode();
    869999        bool change = findLeavingEdge();
     
    8741004          updateDepthPotential();
    8751005        }
    876 #ifdef _DEBUG_ITER_
    877         ++iter_num;
    878 #endif
    879       }
    880 
    881 #ifdef _DEBUG_ITER_
    882       std::cout << "Network Simplex algorithm finished. " << iter_num
    883                 << " pivot iterations performed." << std::endl;
    884 #endif
    885 
    886       // Checking if the flow amount equals zero on all the
    887       // artificial edges
    888       for (InEdgeIt e(graph, root); e != INVALID; ++e)
    889         if (flow[e] > 0) return false;
    890       for (OutEdgeIt e(graph, root); e != INVALID; ++e)
    891         if (flow[e] > 0) return false;
    892 
    893       // Copying flow values to flow_result
    894       if (lower) {
    895         for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
    896           flow_result[e] = (*lower)[e] + flow[edge_ref[e]];
     1006      }
     1007
     1008      // Checking if the flow amount equals zero on all the artificial
     1009      // edges
     1010      for (InEdgeIt e(_graph, _root); e != INVALID; ++e)
     1011        if (_flow[e] > 0) return false;
     1012      for (OutEdgeIt e(_graph, _root); e != INVALID; ++e)
     1013        if (_flow[e] > 0) return false;
     1014
     1015      // Copying flow values to _flow_result
     1016      if (_lower) {
     1017        for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
     1018          _flow_result[e] = (*_lower)[e] + _flow[_edge_ref[e]];
    8971019      } else {
    898         for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
    899           flow_result[e] = flow[edge_ref[e]];
    900       }
    901       // Copying potential values to potential_result
    902       for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
    903         potential_result[n] = potential[node_ref[n]];
     1020        for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
     1021          _flow_result[e] = _flow[_edge_ref[e]];
     1022      }
     1023      // Copying potential values to _potential_result
     1024      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
     1025        _potential_result[n] = _potential[_node_ref[n]];
    9041026
    9051027      return true;
Note: See TracChangeset for help on using the changeset viewer.