lemon/network_simplex.h
author kpeter
Tue, 05 Feb 2008 11:03:55 +0000
changeset 2560 635e7985be46
parent 2553 bfced05fa852
child 2569 12c2c5c4330b
permissions -rw-r--r--
Fix the rename script.
     1 /* -*- C++ -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library
     4  *
     5  * Copyright (C) 2003-2008
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 
    19 #ifndef LEMON_NETWORK_SIMPLEX_H
    20 #define LEMON_NETWORK_SIMPLEX_H
    21 
    22 /// \ingroup min_cost_flow
    23 ///
    24 /// \file
    25 /// \brief The network simplex algorithm for finding a minimum cost flow.
    26 
    27 #include <limits>
    28 #include <lemon/graph_adaptor.h>
    29 #include <lemon/graph_utils.h>
    30 #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 <cmath>
    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
    66 
    67 namespace lemon {
    68 
    69   /// \addtogroup min_cost_flow
    70   /// @{
    71 
    72   /// \brief Implementation of the network simplex algorithm for
    73   /// finding a minimum cost flow.
    74   ///
    75   /// \ref NetworkSimplex implements the network simplex algorithm for
    76   /// finding a minimum cost flow.
    77   ///
    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.
    83   ///
    84   /// \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.
    91   ///
    92   /// \author Peter Kovacs
    93 
    94   template < typename Graph,
    95              typename LowerMap = typename Graph::template EdgeMap<int>,
    96              typename CapacityMap = LowerMap,
    97              typename CostMap = typename Graph::template EdgeMap<int>,
    98              typename SupplyMap = typename Graph::template NodeMap
    99                                   <typename CapacityMap::Value> >
   100   class NetworkSimplex
   101   {
   102     typedef typename LowerMap::Value Lower;
   103     typedef typename CapacityMap::Value Capacity;
   104     typedef typename CostMap::Value Cost;
   105     typedef typename SupplyMap::Value Supply;
   106 
   107     typedef SmartGraph SGraph;
   108     GRAPH_TYPEDEFS(typename SGraph);
   109 
   110     typedef typename SGraph::template EdgeMap<Lower> SLowerMap;
   111     typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap;
   112     typedef typename SGraph::template EdgeMap<Cost> SCostMap;
   113     typedef typename SGraph::template NodeMap<Supply> SSupplyMap;
   114     typedef typename SGraph::template NodeMap<Cost> SPotentialMap;
   115 
   116     typedef typename SGraph::template NodeMap<int> IntNodeMap;
   117     typedef typename SGraph::template NodeMap<bool> BoolNodeMap;
   118     typedef typename SGraph::template NodeMap<Node> NodeNodeMap;
   119     typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap;
   120     typedef typename SGraph::template EdgeMap<int> IntEdgeMap;
   121 
   122     typedef typename Graph::template NodeMap<Node> NodeRefMap;
   123     typedef typename Graph::template EdgeMap<Edge> EdgeRefMap;
   124 
   125   public:
   126 
   127     /// The type of the flow map.
   128     typedef typename Graph::template EdgeMap<Capacity> FlowMap;
   129     /// The type of the potential map.
   130     typedef typename Graph::template NodeMap<Cost> PotentialMap;
   131 
   132   protected:
   133 
   134     /// Map adaptor class for handling reduced edge costs.
   135     class ReducedCostMap : public MapBase<Edge, Cost>
   136     {
   137     private:
   138 
   139       const SGraph &gr;
   140       const SCostMap &cost_map;
   141       const SPotentialMap &pot_map;
   142 
   143     public:
   144 
   145       ReducedCostMap( const SGraph &_gr,
   146                       const SCostMap &_cm,
   147                       const SPotentialMap &_pm ) :
   148         gr(_gr), cost_map(_cm), pot_map(_pm) {}
   149 
   150       Cost operator[](const Edge &e) const {
   151         return cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];
   152       }
   153 
   154     }; //class ReducedCostMap
   155 
   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;
   200 
   201     // The entering edge of the current pivot iteration.
   202     Edge in_edge;
   203     // Temporary nodes used in the current pivot iteration.
   204     Node join, u_in, v_in, u_out, v_out;
   205     Node right, first, second, last;
   206     Node stem, par_stem, new_stem;
   207     // The maximum augment amount along the found cycle in the current
   208     // pivot iteration.
   209     Capacity delta;
   210 
   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 
   238   public :
   239 
   240     /// \brief General constructor of the class (with lower bounds).
   241     ///
   242     /// General constructor of the class (with lower bounds).
   243     ///
   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)
   260     {
   261       // Checking the sum of supply values
   262       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)
   274         .run();
   275 
   276       // 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;
   287       }
   288     }
   289 
   290     /// \brief General constructor of the class (without lower bounds).
   291     ///
   292     /// General constructor of the class (without lower bounds).
   293     ///
   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)
   308     {
   309       // Checking the sum of supply values
   310       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)
   322         .run();
   323     }
   324 
   325     /// \brief Simple constructor of the class (with lower bounds).
   326     ///
   327     /// Simple constructor of the class (with lower bounds).
   328     ///
   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)
   350     {
   351       // Copying graph_ref to graph
   352       copyGraph(graph, graph_ref)
   353         .edgeMap(cost, _cost)
   354         .nodeRef(node_ref)
   355         .edgeRef(edge_ref)
   356         .run();
   357 
   358       // 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;
   373     }
   374 
   375     /// \brief Simple constructor of the class (without lower bounds).
   376     ///
   377     /// Simple constructor of the class (without lower bounds).
   378     ///
   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)
   398     {
   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)
   405         .run();
   406       supply[node_ref[_s]] =  _flow_value;
   407       supply[node_ref[_t]] = -_flow_value;
   408       valid_supply = true;
   409     }
   410 
   411     /// \brief Runs the algorithm.
   412     ///
   413     /// Runs the algorithm.
   414     ///
   415     /// \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.
   423     ///
   424     /// \pre \ref run() must be called before using this function.
   425     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).
   434     ///
   435     /// \pre \ref run() must be called before using this function.
   436     const PotentialMap& potentialMap() const {
   437       return potential_result;
   438     }
   439 
   440     /// \brief Returns the total cost of the found flow.
   441     ///
   442     /// Returns the total cost of the found flow. The complexity of the
   443     /// function is \f$ O(e) \f$.
   444     ///
   445     /// \pre \ref run() must be called before using this function.
   446     Cost totalCost() const {
   447       Cost c = 0;
   448       for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
   449         c += flow_result[e] * cost[edge_ref[e]];
   450       return c;
   451     }
   452 
   453   protected:
   454 
   455     /// \brief Extends the underlaying graph and initializes all the
   456     /// node and edge maps.
   457     bool init() {
   458       if (!valid_supply) return false;
   459 
   460       // Initializing state and flow maps
   461       for (EdgeIt e(graph); e != INVALID; ++e) {
   462         flow[e] = 0;
   463         state[e] = LOWER;
   464       }
   465 
   466       // 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;
   473 
   474       // Adding artificial edges to the graph and initializing the node
   475       // maps of the spanning tree data structure
   476       Supply sum = 0;
   477       Node last = root;
   478       Edge e;
   479       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;
   483         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;
   492         } 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       }
   629       return true;
   630     }
   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.
   679     Node findJoinNode() {
   680       // Finding the join node
   681       Node u = graph.source(in_edge);
   682       Node v = graph.target(in_edge);
   683       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];
   690       }
   691       return u;
   692     }
   693 
   694     /// \brief Finds the leaving edge of the cycle. Returns \c true if
   695     /// the leaving edge is not the same as the entering edge.
   696     bool findLeavingEdge() {
   697       // Initializing first and second nodes according to the direction
   698       // of the cycle
   699       if (state[in_edge] == LOWER) {
   700         first = graph.source(in_edge);
   701         second  = graph.target(in_edge);
   702       } else {
   703         first = graph.target(in_edge);
   704         second  = graph.source(in_edge);
   705       }
   706       delta = capacity[in_edge];
   707       bool result = false;
   708       Capacity d;
   709       Edge e;
   710 
   711       // Searching the cycle along the path form the first node to the
   712       // 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];
   716         if (d < delta) {
   717           delta = d;
   718           u_out = u;
   719           u_in = first;
   720           v_in = second;
   721           result = true;
   722         }
   723       }
   724       // Searching the cycle along the path form the second node to the
   725       // 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];
   729         if (d <= delta) {
   730           delta = d;
   731           u_out = u;
   732           u_in = second;
   733           v_in = first;
   734           result = true;
   735         }
   736       }
   737       return result;
   738     }
   739 
   740     /// \brief Changes \c flow and \c state edge maps.
   741     void changeFlows(bool change) {
   742       // Augmenting along the cycle
   743       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;
   751         }
   752       }
   753       // Updating the state of the entering and leaving edges
   754       if (change) {
   755         state[in_edge] = TREE;
   756         state[pred_edge[u_out]] =
   757           (flow[pred_edge[u_out]] == 0) ? LOWER : UPPER;
   758       } else {
   759         state[in_edge] = -state[in_edge];
   760       }
   761     }
   762 
   763     /// \brief Updates \c thread and \c parent node maps.
   764     void updateThreadParent() {
   765       Node u;
   766       v_out = parent[u_out];
   767 
   768       // Handling the case when join and v_out coincide
   769       bool par_first = false;
   770       if (join == v_out) {
   771         for (u = join; u != u_in && u != v_in; u = thread[u]) ;
   772         if (u == v_in) {
   773           par_first = true;
   774           while (thread[u] != u_out) u = thread[u];
   775           first = u;
   776         }
   777       }
   778 
   779       // Finding the last successor of u_in (u) and the node after it
   780       // (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];
   788 
   789       // Updating stem nodes
   790       thread[v_in] = stem = u_in;
   791       par_stem = v_in;
   792       while (stem != u_out) {
   793         thread[u] = new_stem = parent[stem];
   794 
   795         // Finding the node just before the stem node (u) according to
   796         // the original thread index
   797         for (u = new_stem; thread[u] != stem; u = thread[u]) ;
   798         thread[u] = right;
   799 
   800         // Changing the parent node of stem and shifting stem and
   801         // par_stem nodes
   802         parent[stem] = par_stem;
   803         par_stem = stem;
   804         stem = new_stem;
   805 
   806         // Finding the last successor of stem (u) and the node after it
   807         // (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;
   813 
   814       if (join == v_out && par_first) {
   815         if (first != v_in) thread[first] = right;
   816       } 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.
   823     void updatePredEdge() {
   824       Node u = u_out, v;
   825       while (u != u_in) {
   826         v = parent[u];
   827         pred_edge[u] = pred_edge[v];
   828         forward[u] = !forward[v];
   829         u = v;
   830       }
   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.
   836     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;
   843       while (true) {
   844         v = parent[u];
   845         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.
   856     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       {
   868         join = findJoinNode();
   869         bool change = findLeavingEdge();
   870         changeFlows(change);
   871         if (change) {
   872           updateThreadParent();
   873           updatePredEdge();
   874           updateDepthPotential();
   875         }
   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]];
   897       } 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]];
   904 
   905       return true;
   906     }
   907 
   908   }; //class NetworkSimplex
   909 
   910   ///@}
   911 
   912 } //namespace lemon
   913 
   914 #endif //LEMON_NETWORK_SIMPLEX_H