lemon/network_simplex.h
changeset 2620 8f41a3129746
parent 2593 8eed667ea23c
child 2623 90defb96ee61
equal deleted inserted replaced
13:3dea823e213e 14:104250edf62f
    30 #include <lemon/graph_adaptor.h>
    30 #include <lemon/graph_adaptor.h>
    31 #include <lemon/graph_utils.h>
    31 #include <lemon/graph_utils.h>
    32 #include <lemon/smart_graph.h>
    32 #include <lemon/smart_graph.h>
    33 #include <lemon/math.h>
    33 #include <lemon/math.h>
    34 
    34 
       
    35 #define _DEBUG_
       
    36 
    35 namespace lemon {
    37 namespace lemon {
    36 
    38 
    37   /// \addtogroup min_cost_flow
    39   /// \addtogroup min_cost_flow
    38   /// @{
    40   /// @{
    39 
    41 
    40   /// \brief Implementation of the network simplex algorithm for
    42   /// \brief Implementation of the primal network simplex algorithm
    41   /// finding a minimum cost flow.
    43   /// for finding a minimum cost flow.
    42   ///
    44   ///
    43   /// \ref NetworkSimplex implements the network simplex algorithm for
    45   /// \ref NetworkSimplex implements the primal network simplex algorithm
    44   /// finding a minimum cost flow.
    46   /// for finding a minimum cost flow.
    45   ///
    47   ///
    46   /// \tparam Graph The directed graph type the algorithm runs on.
    48   /// \tparam Graph The directed graph type the algorithm runs on.
    47   /// \tparam LowerMap The type of the lower bound map.
    49   /// \tparam LowerMap The type of the lower bound map.
    48   /// \tparam CapacityMap The type of the capacity (upper bound) map.
    50   /// \tparam CapacityMap The type of the capacity (upper bound) map.
    49   /// \tparam CostMap The type of the cost (length) map.
    51   /// \tparam CostMap The type of the cost (length) map.
    53   /// - Edge capacities and costs should be \e non-negative \e integers.
    55   /// - Edge capacities and costs should be \e non-negative \e integers.
    54   /// - Supply values should be \e signed \e integers.
    56   /// - Supply values should be \e signed \e integers.
    55   /// - The value types of the maps should be convertible to each other.
    57   /// - The value types of the maps should be convertible to each other.
    56   /// - \c CostMap::Value must be signed type.
    58   /// - \c CostMap::Value must be signed type.
    57   ///
    59   ///
    58   /// \note \ref NetworkSimplex provides six different pivot rule
    60   /// \note \ref NetworkSimplex provides five different pivot rule
    59   /// implementations that significantly affect the efficiency of the
    61   /// implementations that significantly affect the efficiency of the
    60   /// algorithm.
    62   /// algorithm.
    61   /// By default a combined pivot rule is used, which is the fastest
    63   /// By default "Block Search" pivot rule is used, which proved to be
    62   /// implementation according to our benchmark tests.
    64   /// by far the most efficient according to our benchmark tests.
    63   /// Another pivot rule can be selected using \ref run() function
    65   /// However another pivot rule can be selected using \ref run()
    64   /// with the proper parameter.
    66   /// function with the proper parameter.
    65   ///
    67   ///
    66   /// \author Peter Kovacs
    68   /// \author Peter Kovacs
    67 
       
    68   template < typename Graph,
    69   template < typename Graph,
    69              typename LowerMap = typename Graph::template EdgeMap<int>,
    70              typename LowerMap = typename Graph::template EdgeMap<int>,
    70              typename CapacityMap = typename Graph::template EdgeMap<int>,
    71              typename CapacityMap = typename Graph::template EdgeMap<int>,
    71              typename CostMap = typename Graph::template EdgeMap<int>,
    72              typename CostMap = typename Graph::template EdgeMap<int>,
    72              typename SupplyMap = typename Graph::template NodeMap<int> >
    73              typename SupplyMap = typename Graph::template NodeMap<int> >
    87     typedef typename SGraph::template NodeMap<int> IntNodeMap;
    88     typedef typename SGraph::template NodeMap<int> IntNodeMap;
    88     typedef typename SGraph::template NodeMap<bool> BoolNodeMap;
    89     typedef typename SGraph::template NodeMap<bool> BoolNodeMap;
    89     typedef typename SGraph::template NodeMap<Node> NodeNodeMap;
    90     typedef typename SGraph::template NodeMap<Node> NodeNodeMap;
    90     typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap;
    91     typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap;
    91     typedef typename SGraph::template EdgeMap<int> IntEdgeMap;
    92     typedef typename SGraph::template EdgeMap<int> IntEdgeMap;
       
    93     typedef typename SGraph::template EdgeMap<bool> BoolEdgeMap;
    92 
    94 
    93     typedef typename Graph::template NodeMap<Node> NodeRefMap;
    95     typedef typename Graph::template NodeMap<Node> NodeRefMap;
    94     typedef typename Graph::template EdgeMap<Edge> EdgeRefMap;
    96     typedef typename Graph::template EdgeMap<Edge> EdgeRefMap;
       
    97 
       
    98     typedef std::vector<Edge> EdgeVector;
    95 
    99 
    96   public:
   100   public:
    97 
   101 
    98     /// The type of the flow map.
   102     /// The type of the flow map.
    99     typedef typename Graph::template EdgeMap<Capacity> FlowMap;
   103     typedef typename Graph::template EdgeMap<Capacity> FlowMap;
   105     /// Enum type to select the pivot rule used by \ref run().
   109     /// Enum type to select the pivot rule used by \ref run().
   106     enum PivotRuleEnum {
   110     enum PivotRuleEnum {
   107       FIRST_ELIGIBLE_PIVOT,
   111       FIRST_ELIGIBLE_PIVOT,
   108       BEST_ELIGIBLE_PIVOT,
   112       BEST_ELIGIBLE_PIVOT,
   109       BLOCK_SEARCH_PIVOT,
   113       BLOCK_SEARCH_PIVOT,
   110       LIMITED_SEARCH_PIVOT,
       
   111       CANDIDATE_LIST_PIVOT,
   114       CANDIDATE_LIST_PIVOT,
   112       COMBINED_PIVOT
   115       ALTERING_LIST_PIVOT
   113     };
   116     };
   114 
   117 
   115   private:
   118   private:
   116 
   119 
   117     /// \brief Map adaptor class for handling reduced edge costs.
   120     /// \brief Map adaptor class for handling reduced edge costs.
   146     /// \brief Implementation of the "First Eligible" pivot rule for the
   149     /// \brief Implementation of the "First Eligible" pivot rule for the
   147     /// \ref NetworkSimplex "network simplex" algorithm.
   150     /// \ref NetworkSimplex "network simplex" algorithm.
   148     ///
   151     ///
   149     /// This class implements the "First Eligible" pivot rule
   152     /// This class implements the "First Eligible" pivot rule
   150     /// for the \ref NetworkSimplex "network simplex" algorithm.
   153     /// for the \ref NetworkSimplex "network simplex" algorithm.
       
   154     ///
       
   155     /// For more information see \ref NetworkSimplex::run().
   151     class FirstEligiblePivotRule
   156     class FirstEligiblePivotRule
   152     {
   157     {
   153     private:
   158     private:
   154 
   159 
       
   160       // References to the NetworkSimplex class
   155       NetworkSimplex &_ns;
   161       NetworkSimplex &_ns;
   156       EdgeIt _next_edge;
   162       EdgeVector &_edges;
       
   163 
       
   164       int _next_edge;
   157 
   165 
   158     public:
   166     public:
   159 
   167 
   160       /// Constructor.
   168       /// Constructor
   161       FirstEligiblePivotRule(NetworkSimplex &ns) :
   169       FirstEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) :
   162         _ns(ns), _next_edge(ns._graph) {}
   170         _ns(ns), _edges(edges), _next_edge(0) {}
   163 
   171 
   164       /// Finds the next entering edge.
   172       /// Find next entering edge
   165       bool findEnteringEdge() {
   173       inline bool findEnteringEdge() {
   166         for (EdgeIt e = _next_edge; e != INVALID; ++e) {
   174         Edge e;
       
   175         for (int i = _next_edge; i < int(_edges.size()); ++i) {
       
   176           e = _edges[i];
   167           if (_ns._state[e] * _ns._red_cost[e] < 0) {
   177           if (_ns._state[e] * _ns._red_cost[e] < 0) {
   168             _ns._in_edge = e;
   178             _ns._in_edge = e;
   169             _next_edge = ++e;
   179             _next_edge = i + 1;
   170             return true;
   180             return true;
   171           }
   181           }
   172         }
   182         }
   173         for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
   183         for (int i = 0; i < _next_edge; ++i) {
       
   184           e = _edges[i];
   174           if (_ns._state[e] * _ns._red_cost[e] < 0) {
   185           if (_ns._state[e] * _ns._red_cost[e] < 0) {
   175             _ns._in_edge = e;
   186             _ns._in_edge = e;
   176             _next_edge = ++e;
   187             _next_edge = i + 1;
   177             return true;
   188             return true;
   178           }
   189           }
   179         }
   190         }
   180         return false;
   191         return false;
   181       }
   192       }
   184     /// \brief Implementation of the "Best Eligible" pivot rule for the
   195     /// \brief Implementation of the "Best Eligible" pivot rule for the
   185     /// \ref NetworkSimplex "network simplex" algorithm.
   196     /// \ref NetworkSimplex "network simplex" algorithm.
   186     ///
   197     ///
   187     /// This class implements the "Best Eligible" pivot rule
   198     /// This class implements the "Best Eligible" pivot rule
   188     /// for the \ref NetworkSimplex "network simplex" algorithm.
   199     /// for the \ref NetworkSimplex "network simplex" algorithm.
       
   200     ///
       
   201     /// For more information see \ref NetworkSimplex::run().
   189     class BestEligiblePivotRule
   202     class BestEligiblePivotRule
   190     {
   203     {
   191     private:
   204     private:
   192 
   205 
       
   206       // References to the NetworkSimplex class
   193       NetworkSimplex &_ns;
   207       NetworkSimplex &_ns;
       
   208       EdgeVector &_edges;
   194 
   209 
   195     public:
   210     public:
   196 
   211 
   197       /// Constructor.
   212       /// Constructor
   198       BestEligiblePivotRule(NetworkSimplex &ns) : _ns(ns) {}
   213       BestEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) :
   199 
   214         _ns(ns), _edges(edges) {}
   200       /// Finds the next entering edge.
   215 
   201       bool findEnteringEdge() {
   216       /// Find next entering edge
       
   217       inline bool findEnteringEdge() {
   202         Cost min = 0;
   218         Cost min = 0;
   203         for (EdgeIt e(_ns._graph); e != INVALID; ++e) {
   219         Edge e;
       
   220         for (int i = 0; i < int(_edges.size()); ++i) {
       
   221           e = _edges[i];
   204           if (_ns._state[e] * _ns._red_cost[e] < min) {
   222           if (_ns._state[e] * _ns._red_cost[e] < min) {
   205             min = _ns._state[e] * _ns._red_cost[e];
   223             min = _ns._state[e] * _ns._red_cost[e];
   206             _ns._in_edge = e;
   224             _ns._in_edge = e;
   207           }
   225           }
   208         }
   226         }
   213     /// \brief Implementation of the "Block Search" pivot rule for the
   231     /// \brief Implementation of the "Block Search" pivot rule for the
   214     /// \ref NetworkSimplex "network simplex" algorithm.
   232     /// \ref NetworkSimplex "network simplex" algorithm.
   215     ///
   233     ///
   216     /// This class implements the "Block Search" pivot rule
   234     /// This class implements the "Block Search" pivot rule
   217     /// for the \ref NetworkSimplex "network simplex" algorithm.
   235     /// for the \ref NetworkSimplex "network simplex" algorithm.
       
   236     ///
       
   237     /// For more information see \ref NetworkSimplex::run().
   218     class BlockSearchPivotRule
   238     class BlockSearchPivotRule
   219     {
   239     {
   220     private:
   240     private:
   221 
   241 
       
   242       // References to the NetworkSimplex class
   222       NetworkSimplex &_ns;
   243       NetworkSimplex &_ns;
   223       EdgeIt _next_edge, _min_edge;
   244       EdgeVector &_edges;
       
   245 
   224       int _block_size;
   246       int _block_size;
   225 
   247       int _next_edge, _min_edge;
   226       static const int MIN_BLOCK_SIZE = 10;
       
   227 
   248 
   228     public:
   249     public:
   229 
   250 
   230       /// Constructor.
   251       /// Constructor
   231       BlockSearchPivotRule(NetworkSimplex &ns) :
   252       BlockSearchPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
   232         _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
   253         _ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
   233       {
   254       {
   234         _block_size = 2 * int(sqrt(countEdges(_ns._graph)));
   255         // The main parameters of the pivot rule
   235         if (_block_size < MIN_BLOCK_SIZE) _block_size = MIN_BLOCK_SIZE;
   256         const double BLOCK_SIZE_FACTOR = 2.0;
   236       }
   257         const int MIN_BLOCK_SIZE = 10;
   237 
   258 
   238       /// Finds the next entering edge.
   259         _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
   239       bool findEnteringEdge() {
   260                                 MIN_BLOCK_SIZE );
       
   261       }
       
   262 
       
   263       /// Find next entering edge
       
   264       inline bool findEnteringEdge() {
   240         Cost curr, min = 0;
   265         Cost curr, min = 0;
   241         int cnt = 0;
   266         Edge e;
   242         for (EdgeIt e = _next_edge; e != INVALID; ++e) {
   267         int cnt = _block_size;
       
   268         int i;
       
   269         for (i = _next_edge; i < int(_edges.size()); ++i) {
       
   270           e = _edges[i];
   243           if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
   271           if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
   244             min = curr;
   272             min = curr;
   245             _min_edge = e;
   273             _min_edge = i;
   246           }
   274           }
   247           if (++cnt == _block_size) {
   275           if (--cnt == 0) {
   248             if (min < 0) break;
   276             if (min < 0) break;
   249             cnt = 0;
   277             cnt = _block_size;
   250           }
   278           }
   251         }
   279         }
   252         if (min == 0) {
   280         if (min == 0 || cnt > 0) {
   253           for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
   281           for (i = 0; i < _next_edge; ++i) {
       
   282             e = _edges[i];
   254             if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
   283             if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
   255               min = curr;
   284               min = curr;
   256               _min_edge = e;
   285               _min_edge = i;
   257             }
   286             }
   258             if (++cnt == _block_size) {
   287             if (--cnt == 0) {
   259               if (min < 0) break;
   288               if (min < 0) break;
   260               cnt = 0;
   289               cnt = _block_size;
   261             }
   290             }
   262           }
   291           }
   263         }
   292         }
   264         _ns._in_edge = _min_edge;
   293         if (min >= 0) return false;
   265         _next_edge = ++_min_edge;
   294         _ns._in_edge = _edges[_min_edge];
   266         return min < 0;
   295         _next_edge = i;
       
   296         return true;
   267       }
   297       }
   268     }; //class BlockSearchPivotRule
   298     }; //class BlockSearchPivotRule
   269 
       
   270     /// \brief Implementation of the "Limited Search" pivot rule for the
       
   271     /// \ref NetworkSimplex "network simplex" algorithm.
       
   272     ///
       
   273     /// This class implements the "Limited Search" pivot rule
       
   274     /// for the \ref NetworkSimplex "network simplex" algorithm.
       
   275     class LimitedSearchPivotRule
       
   276     {
       
   277     private:
       
   278 
       
   279       NetworkSimplex &_ns;
       
   280       EdgeIt _next_edge, _min_edge;
       
   281       int _sample_size;
       
   282 
       
   283       static const int SAMPLE_SIZE_FACTOR = 15;
       
   284       static const int MIN_SAMPLE_SIZE = 10;
       
   285 
       
   286     public:
       
   287 
       
   288       /// Constructor.
       
   289       LimitedSearchPivotRule(NetworkSimplex &ns) :
       
   290         _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
       
   291       {
       
   292         _sample_size = countEdges(_ns._graph) *
       
   293                        SAMPLE_SIZE_FACTOR / 10000;
       
   294         if (_sample_size < MIN_SAMPLE_SIZE)
       
   295           _sample_size = MIN_SAMPLE_SIZE;
       
   296       }
       
   297 
       
   298       /// Finds the next entering edge.
       
   299       bool findEnteringEdge() {
       
   300         Cost curr, min = 0;
       
   301         int cnt = 0;
       
   302         for (EdgeIt e = _next_edge; e != INVALID; ++e) {
       
   303           if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
       
   304             min = curr;
       
   305             _min_edge = e;
       
   306           }
       
   307           if (curr < 0 && ++cnt == _sample_size) break;
       
   308         }
       
   309         if (min == 0) {
       
   310           for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
       
   311             if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
       
   312               min = curr;
       
   313               _min_edge = e;
       
   314             }
       
   315             if (curr < 0 && ++cnt == _sample_size) break;
       
   316           }
       
   317         }
       
   318         _ns._in_edge = _min_edge;
       
   319         _next_edge = ++_min_edge;
       
   320         return min < 0;
       
   321       }
       
   322     }; //class LimitedSearchPivotRule
       
   323 
   299 
   324     /// \brief Implementation of the "Candidate List" pivot rule for the
   300     /// \brief Implementation of the "Candidate List" pivot rule for the
   325     /// \ref NetworkSimplex "network simplex" algorithm.
   301     /// \ref NetworkSimplex "network simplex" algorithm.
   326     ///
   302     ///
   327     /// This class implements the "Candidate List" pivot rule
   303     /// This class implements the "Candidate List" pivot rule
   328     /// for the \ref NetworkSimplex "network simplex" algorithm.
   304     /// for the \ref NetworkSimplex "network simplex" algorithm.
       
   305     ///
       
   306     /// For more information see \ref NetworkSimplex::run().
   329     class CandidateListPivotRule
   307     class CandidateListPivotRule
   330     {
   308     {
   331     private:
   309     private:
   332 
   310 
       
   311       // References to the NetworkSimplex class
   333       NetworkSimplex &_ns;
   312       NetworkSimplex &_ns;
   334 
   313       EdgeVector &_edges;
   335       // The list of candidate edges.
   314 
   336       std::vector<Edge> _candidates;
   315       EdgeVector _candidates;
   337       // The maximum length of the edge list.
   316       int _list_length, _minor_limit;
   338       int _list_length;
   317       int _curr_length, _minor_count;
   339       // The maximum number of minor iterations between two major
   318       int _next_edge, _min_edge;
   340       // itarations.
       
   341       int _minor_limit;
       
   342 
       
   343       int _minor_count;
       
   344       EdgeIt _next_edge;
       
   345 
       
   346       static const int LIST_LENGTH_FACTOR = 20;
       
   347       static const int MINOR_LIMIT_FACTOR = 10;
       
   348       static const int MIN_LIST_LENGTH = 10;
       
   349       static const int MIN_MINOR_LIMIT = 2;
       
   350 
   319 
   351     public:
   320     public:
   352 
   321 
   353       /// Constructor.
   322       /// Constructor
   354       CandidateListPivotRule(NetworkSimplex &ns) :
   323       CandidateListPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
   355         _ns(ns), _next_edge(ns._graph)
   324         _ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
   356       {
   325       {
   357         int edge_num = countEdges(_ns._graph);
   326         // The main parameters of the pivot rule
   358         _minor_count = 0;
   327         const double LIST_LENGTH_FACTOR = 1.0;
   359         _list_length = edge_num * LIST_LENGTH_FACTOR / 10000;
   328         const int MIN_LIST_LENGTH = 10;
   360         if (_list_length < MIN_LIST_LENGTH)
   329         const double MINOR_LIMIT_FACTOR = 0.1;
   361           _list_length = MIN_LIST_LENGTH;
   330         const int MIN_MINOR_LIMIT = 3;
   362         _minor_limit = _list_length * MINOR_LIMIT_FACTOR / 100;
   331 
   363         if (_minor_limit < MIN_MINOR_LIMIT)
   332         _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edges.size())),
   364           _minor_limit = MIN_MINOR_LIMIT;
   333                                  MIN_LIST_LENGTH );
   365       }
   334         _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
   366 
   335                                  MIN_MINOR_LIMIT );
   367       /// Finds the next entering edge.
   336         _curr_length = _minor_count = 0;
   368       bool findEnteringEdge() {
   337         _candidates.resize(_list_length);
       
   338       }
       
   339 
       
   340       /// Find next entering edge
       
   341       inline bool findEnteringEdge() {
   369         Cost min, curr;
   342         Cost min, curr;
   370         if (_minor_count < _minor_limit && _candidates.size() > 0) {
   343         if (_curr_length > 0 && _minor_count < _minor_limit) {
   371           // Minor iteration
   344           // Minor iteration: selecting the best eligible edge from
       
   345           // the current candidate list
   372           ++_minor_count;
   346           ++_minor_count;
   373           Edge e;
   347           Edge e;
   374           min = 0;
   348           min = 0;
   375           for (int i = 0; i < int(_candidates.size()); ++i) {
   349           for (int i = 0; i < _curr_length; ++i) {
   376             e = _candidates[i];
   350             e = _candidates[i];
   377             if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
   351             curr = _ns._state[e] * _ns._red_cost[e];
   378               min = curr;
       
   379               _ns._in_edge = e;
       
   380             }
       
   381           }
       
   382           if (min < 0) return true;
       
   383         }
       
   384 
       
   385         // Major iteration
       
   386         _candidates.clear();
       
   387         EdgeIt e = _next_edge;
       
   388         min = 0;
       
   389         for ( ; e != INVALID; ++e) {
       
   390           if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
       
   391             _candidates.push_back(e);
       
   392             if (curr < min) {
   352             if (curr < min) {
   393               min = curr;
   353               min = curr;
   394               _ns._in_edge = e;
   354               _ns._in_edge = e;
   395             }
   355             }
   396             if (int(_candidates.size()) == _list_length) break;
   356             if (curr >= 0) {
       
   357               _candidates[i--] = _candidates[--_curr_length];
       
   358             }
   397           }
   359           }
   398         }
   360           if (min < 0) return true;
   399         if (int(_candidates.size()) < _list_length) {
   361         }
   400           for (e = EdgeIt(_ns._graph); e != _next_edge; ++e) {
   362 
       
   363         // Major iteration: building a new candidate list
       
   364         Edge e;
       
   365         min = 0;
       
   366         _curr_length = 0;
       
   367         int i;
       
   368         for (i = _next_edge; i < int(_edges.size()); ++i) {
       
   369           e = _edges[i];
       
   370           if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
       
   371             _candidates[_curr_length++] = e;
       
   372             if (curr < min) {
       
   373               min = curr;
       
   374               _min_edge = i;
       
   375             }
       
   376             if (_curr_length == _list_length) break;
       
   377           }
       
   378         }
       
   379         if (_curr_length < _list_length) {
       
   380           for (i = 0; i < _next_edge; ++i) {
       
   381             e = _edges[i];
   401             if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
   382             if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
   402               _candidates.push_back(e);
   383               _candidates[_curr_length++] = e;
   403               if (curr < min) {
   384               if (curr < min) {
   404                 min = curr;
   385                 min = curr;
   405                 _ns._in_edge = e;
   386                 _min_edge = i;
   406               }
   387               }
   407               if (int(_candidates.size()) == _list_length) break;
   388               if (_curr_length == _list_length) break;
   408             }
   389             }
   409           }
   390           }
   410         }
   391         }
   411         if (_candidates.size() == 0) return false;
   392         if (_curr_length == 0) return false;
   412         _minor_count = 1;
   393         _minor_count = 1;
   413         _next_edge = ++e;
   394         _ns._in_edge = _edges[_min_edge];
       
   395         _next_edge = i;
   414         return true;
   396         return true;
   415       }
   397       }
   416     }; //class CandidateListPivotRule
   398     }; //class CandidateListPivotRule
       
   399 
       
   400     /// \brief Implementation of the "Altering Candidate List" pivot rule
       
   401     /// for the \ref NetworkSimplex "network simplex" algorithm.
       
   402     ///
       
   403     /// This class implements the "Altering Candidate List" pivot rule
       
   404     /// for the \ref NetworkSimplex "network simplex" algorithm.
       
   405     ///
       
   406     /// For more information see \ref NetworkSimplex::run().
       
   407     class AlteringListPivotRule
       
   408     {
       
   409     private:
       
   410 
       
   411       // References to the NetworkSimplex class
       
   412       NetworkSimplex &_ns;
       
   413       EdgeVector &_edges;
       
   414 
       
   415       EdgeVector _candidates;
       
   416       SCostMap _cand_cost;
       
   417       int _block_size, _head_length, _curr_length;
       
   418       int _next_edge;
       
   419 
       
   420       // Functor class to compare edges during sort of the candidate list
       
   421       class SortFunc
       
   422       {
       
   423       private:
       
   424         const SCostMap &_map;
       
   425       public:
       
   426         SortFunc(const SCostMap &map) : _map(map) {}
       
   427         bool operator()(const Edge &e1, const Edge &e2) {
       
   428 	  return _map[e1] < _map[e2];
       
   429         }
       
   430       };
       
   431 
       
   432       SortFunc _sort_func;
       
   433 
       
   434     public:
       
   435 
       
   436       /// Constructor
       
   437       AlteringListPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
       
   438         _ns(ns), _edges(edges), _cand_cost(_ns._graph),
       
   439         _next_edge(0), _sort_func(_cand_cost)
       
   440       {
       
   441         // The main parameters of the pivot rule
       
   442         const double BLOCK_SIZE_FACTOR = 1.0;
       
   443         const int MIN_BLOCK_SIZE = 10;
       
   444         const double HEAD_LENGTH_FACTOR = 0.1;
       
   445         const int MIN_HEAD_LENGTH = 5;
       
   446 
       
   447         _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
       
   448                                 MIN_BLOCK_SIZE );
       
   449         _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
       
   450                                  MIN_HEAD_LENGTH );
       
   451         _candidates.resize(_head_length + _block_size);
       
   452         _curr_length = 0;
       
   453       }
       
   454 
       
   455       /// Find next entering edge
       
   456       inline bool findEnteringEdge() {
       
   457         // Checking the current candidate list
       
   458         Edge e;
       
   459         for (int idx = 0; idx < _curr_length; ++idx) {
       
   460           e = _candidates[idx];
       
   461           if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) >= 0) {
       
   462             _candidates[idx--] = _candidates[--_curr_length];
       
   463           }
       
   464         }
       
   465 
       
   466         // Extending the list
       
   467         int cnt = _block_size;
       
   468         int last_edge = 0;
       
   469         int limit = _head_length;
       
   470         for (int i = _next_edge; i < int(_edges.size()); ++i) {
       
   471           e = _edges[i];
       
   472           if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
       
   473             _candidates[_curr_length++] = e;
       
   474             last_edge = i;
       
   475           }
       
   476           if (--cnt == 0) {
       
   477             if (_curr_length > limit) break;
       
   478             limit = 0;
       
   479             cnt = _block_size;
       
   480           }
       
   481         }
       
   482         if (_curr_length <= limit) {
       
   483           for (int i = 0; i < _next_edge; ++i) {
       
   484             e = _edges[i];
       
   485             if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
       
   486               _candidates[_curr_length++] = e;
       
   487               last_edge = i;
       
   488             }
       
   489             if (--cnt == 0) {
       
   490               if (_curr_length > limit) break;
       
   491               limit = 0;
       
   492               cnt = _block_size;
       
   493             }
       
   494           }
       
   495         }
       
   496         if (_curr_length == 0) return false;
       
   497         _next_edge = last_edge + 1;
       
   498 
       
   499         // Sorting the list partially
       
   500         EdgeVector::iterator sort_end = _candidates.begin();
       
   501         EdgeVector::iterator vector_end = _candidates.begin();
       
   502         for (int idx = 0; idx < _curr_length; ++idx) {
       
   503           ++vector_end;
       
   504           if (idx <= _head_length) ++sort_end;
       
   505         }
       
   506         partial_sort(_candidates.begin(), sort_end, vector_end, _sort_func);
       
   507 
       
   508         _ns._in_edge = _candidates[0];
       
   509         if (_curr_length > _head_length) {
       
   510           _candidates[0] = _candidates[_head_length - 1];
       
   511           _curr_length = _head_length - 1;
       
   512         } else {
       
   513           _candidates[0] = _candidates[_curr_length - 1];
       
   514           --_curr_length;
       
   515         }
       
   516 
       
   517         return true;
       
   518       }
       
   519     }; //class AlteringListPivotRule
   417 
   520 
   418   private:
   521   private:
   419 
   522 
   420     // State constants for edges
   523     // State constants for edges
   421     enum EdgeStateEnum {
   524     enum EdgeStateEnum {
   422       STATE_UPPER = -1,
   525       STATE_UPPER = -1,
   423       STATE_TREE  =  0,
   526       STATE_TREE  =  0,
   424       STATE_LOWER =  1
   527       STATE_LOWER =  1
   425     };
   528     };
   426 
       
   427     // Constant for the combined pivot rule.
       
   428     static const int COMBINED_PIVOT_MAX_DEG = 5;
       
   429 
   529 
   430   private:
   530   private:
   431 
   531 
   432     // The directed graph the algorithm runs on
   532     // The directed graph the algorithm runs on
   433     SGraph _graph;
   533     SGraph _graph;
   463     // The root node of the starting spanning tree
   563     // The root node of the starting spanning tree
   464     Node _root;
   564     Node _root;
   465 
   565 
   466     // The reduced cost map
   566     // The reduced cost map
   467     ReducedCostMap _red_cost;
   567     ReducedCostMap _red_cost;
       
   568 
       
   569     // The non-artifical edges
       
   570     EdgeVector _edges;
   468 
   571 
   469     // Members for handling the original graph
   572     // Members for handling the original graph
   470     FlowMap *_flow_result;
   573     FlowMap *_flow_result;
   471     PotentialMap *_potential_result;
   574     PotentialMap *_potential_result;
   472     bool _local_flow;
   575     bool _local_flow;
   670     ~NetworkSimplex() {
   773     ~NetworkSimplex() {
   671       if (_local_flow) delete _flow_result;
   774       if (_local_flow) delete _flow_result;
   672       if (_local_potential) delete _potential_result;
   775       if (_local_potential) delete _potential_result;
   673     }
   776     }
   674 
   777 
   675     /// \brief Sets the flow map.
   778     /// \brief Set the flow map.
   676     ///
   779     ///
   677     /// Sets the flow map.
   780     /// Set the flow map.
   678     ///
   781     ///
   679     /// \return \c (*this)
   782     /// \return \c (*this)
   680     NetworkSimplex& flowMap(FlowMap &map) {
   783     NetworkSimplex& flowMap(FlowMap &map) {
   681       if (_local_flow) {
   784       if (_local_flow) {
   682         delete _flow_result;
   785         delete _flow_result;
   684       }
   787       }
   685       _flow_result = &map;
   788       _flow_result = &map;
   686       return *this;
   789       return *this;
   687     }
   790     }
   688 
   791 
   689     /// \brief Sets the potential map.
   792     /// \brief Set the potential map.
   690     ///
   793     ///
   691     /// Sets the potential map.
   794     /// Set the potential map.
   692     ///
   795     ///
   693     /// \return \c (*this)
   796     /// \return \c (*this)
   694     NetworkSimplex& potentialMap(PotentialMap &map) {
   797     NetworkSimplex& potentialMap(PotentialMap &map) {
   695       if (_local_potential) {
   798       if (_local_potential) {
   696         delete _potential_result;
   799         delete _potential_result;
   699       _potential_result = &map;
   802       _potential_result = &map;
   700       return *this;
   803       return *this;
   701     }
   804     }
   702 
   805 
   703     /// \name Execution control
   806     /// \name Execution control
   704     /// The only way to execute the algorithm is to call the run()
       
   705     /// function.
       
   706 
   807 
   707     /// @{
   808     /// @{
   708 
   809 
   709     /// \brief Runs the algorithm.
   810     /// \brief Runs the algorithm.
   710     ///
   811     ///
   724     ///
   825     ///
   725     /// - BLOCK_SEARCH_PIVOT A specified number of edges are examined in
   826     /// - BLOCK_SEARCH_PIVOT A specified number of edges are examined in
   726     /// every iteration in a wraparound fashion and the best eligible
   827     /// every iteration in a wraparound fashion and the best eligible
   727     /// edge is selected from this block (\ref BlockSearchPivotRule).
   828     /// edge is selected from this block (\ref BlockSearchPivotRule).
   728     ///
   829     ///
   729     /// - LIMITED_SEARCH_PIVOT A specified number of eligible edges are
   830     /// - CANDIDATE_LIST_PIVOT In a major iteration a candidate list is
   730     /// examined in every iteration in a wraparound fashion and the best
   831     /// built from eligible edges in a wraparound fashion and in the
   731     /// one is selected from them (\ref LimitedSearchPivotRule).
   832     /// following minor iterations the best eligible edge is selected
   732     ///
   833     /// from this list (\ref CandidateListPivotRule).
   733     /// - CANDIDATE_LIST_PIVOT In major iterations a candidate list is
   834     ///
   734     /// built from eligible edges and it is used for edge selection in
   835     /// - ALTERING_LIST_PIVOT It is a modified version of the
   735     /// the following minor iterations (\ref CandidateListPivotRule).
   836     /// "Candidate List" pivot rule. It keeps only the several best
   736     ///
   837     /// eligible edges from the former candidate list and extends this
   737     /// - COMBINED_PIVOT This is a combined version of the two fastest
   838     /// list in every iteration (\ref AlteringListPivotRule).
   738     /// pivot rules.
   839     ///
   739     /// For rather sparse graphs \ref LimitedSearchPivotRule
   840     /// According to our comprehensive benchmark tests the "Block Search"
   740     /// "Limited Search" implementation is used, otherwise
   841     /// pivot rule proved to be by far the fastest and the most robust
   741     /// \ref BlockSearchPivotRule "Block Search" pivot rule is used.
   842     /// on various test inputs. Thus it is the default option.
   742     /// According to our benchmark tests this combined method is the
       
   743     /// most efficient.
       
   744     ///
   843     ///
   745     /// \return \c true if a feasible flow can be found.
   844     /// \return \c true if a feasible flow can be found.
   746     bool run(PivotRuleEnum pivot_rule = COMBINED_PIVOT) {
   845     bool run(PivotRuleEnum pivot_rule = BLOCK_SEARCH_PIVOT) {
   747       return init() && start(pivot_rule);
   846       return init() && start(pivot_rule);
   748     }
   847     }
   749 
   848 
   750     /// @}
   849     /// @}
   751 
   850 
   752     /// \name Query Functions
   851     /// \name Query Functions
   753     /// The result of the algorithm can be obtained using these
   852     /// The results of the algorithm can be obtained using these
   754     /// functions.
   853     /// functions.\n
   755     /// \n run() must be called before using them.
   854     /// \ref lemon::NetworkSimplex::run() "run()" must be called before
       
   855     /// using them.
   756 
   856 
   757     /// @{
   857     /// @{
   758 
   858 
   759     /// \brief Returns a const reference to the edge map storing the
   859     /// \brief Return a const reference to the edge map storing the
   760     /// found flow.
   860     /// found flow.
   761     ///
   861     ///
   762     /// Returns a const reference to the edge map storing the found flow.
   862     /// Return a const reference to the edge map storing the found flow.
   763     ///
   863     ///
   764     /// \pre \ref run() must be called before using this function.
   864     /// \pre \ref run() must be called before using this function.
   765     const FlowMap& flowMap() const {
   865     const FlowMap& flowMap() const {
   766       return *_flow_result;
   866       return *_flow_result;
   767     }
   867     }
   768 
   868 
   769     /// \brief Returns a const reference to the node map storing the
   869     /// \brief Return a const reference to the node map storing the
   770     /// found potentials (the dual solution).
   870     /// found potentials (the dual solution).
   771     ///
   871     ///
   772     /// Returns a const reference to the node map storing the found
   872     /// Return a const reference to the node map storing the found
   773     /// potentials (the dual solution).
   873     /// potentials (the dual solution).
   774     ///
   874     ///
   775     /// \pre \ref run() must be called before using this function.
   875     /// \pre \ref run() must be called before using this function.
   776     const PotentialMap& potentialMap() const {
   876     const PotentialMap& potentialMap() const {
   777       return *_potential_result;
   877       return *_potential_result;
   778     }
   878     }
   779 
   879 
   780     /// \brief Returns the flow on the given edge.
   880     /// \brief Return the flow on the given edge.
   781     ///
   881     ///
   782     /// Returns the flow on the given edge.
   882     /// Return the flow on the given edge.
   783     ///
   883     ///
   784     /// \pre \ref run() must be called before using this function.
   884     /// \pre \ref run() must be called before using this function.
   785     Capacity flow(const typename Graph::Edge& edge) const {
   885     Capacity flow(const typename Graph::Edge& edge) const {
   786       return (*_flow_result)[edge];
   886       return (*_flow_result)[edge];
   787     }
   887     }
   788 
   888 
   789     /// \brief Returns the potential of the given node.
   889     /// \brief Return the potential of the given node.
   790     ///
   890     ///
   791     /// Returns the potential of the given node.
   891     /// Return the potential of the given node.
   792     ///
   892     ///
   793     /// \pre \ref run() must be called before using this function.
   893     /// \pre \ref run() must be called before using this function.
   794     Cost potential(const typename Graph::Node& node) const {
   894     Cost potential(const typename Graph::Node& node) const {
   795       return (*_potential_result)[node];
   895       return (*_potential_result)[node];
   796     }
   896     }
   797 
   897 
   798     /// \brief Returns the total cost of the found flow.
   898     /// \brief Return the total cost of the found flow.
   799     ///
   899     ///
   800     /// Returns the total cost of the found flow. The complexity of the
   900     /// Return the total cost of the found flow. The complexity of the
   801     /// function is \f$ O(e) \f$.
   901     /// function is \f$ O(e) \f$.
   802     ///
   902     ///
   803     /// \pre \ref run() must be called before using this function.
   903     /// \pre \ref run() must be called before using this function.
   804     Cost totalCost() const {
   904     Cost totalCost() const {
   805       Cost c = 0;
   905       Cost c = 0;
   810 
   910 
   811     /// @}
   911     /// @}
   812 
   912 
   813   private:
   913   private:
   814 
   914 
   815     /// \brief Extends the underlaying graph and initializes all the
   915     /// \brief Extend the underlying graph and initialize all the
   816     /// node and edge maps.
   916     /// node and edge maps.
   817     bool init() {
   917     bool init() {
   818       if (!_valid_supply) return false;
   918       if (!_valid_supply) return false;
   819 
   919 
   820       // Initializing result maps
   920       // Initializing result maps
   825       if (!_potential_result) {
   925       if (!_potential_result) {
   826         _potential_result = new PotentialMap(_graph_ref);
   926         _potential_result = new PotentialMap(_graph_ref);
   827         _local_potential = true;
   927         _local_potential = true;
   828       }
   928       }
   829 
   929 
   830       // Initializing state and flow maps
   930       // Initializing the edge vector and the edge maps
       
   931       _edges.reserve(countEdges(_graph));
   831       for (EdgeIt e(_graph); e != INVALID; ++e) {
   932       for (EdgeIt e(_graph); e != INVALID; ++e) {
       
   933         _edges.push_back(e);
   832         _flow[e] = 0;
   934         _flow[e] = 0;
   833         _state[e] = STATE_LOWER;
   935         _state[e] = STATE_LOWER;
   834       }
   936       }
   835 
   937 
   836       // Adding an artificial root node to the graph
   938       // Adding an artificial root node to the graph
   871       _thread[last] = _root;
   973       _thread[last] = _root;
   872 
   974 
   873       return true;
   975       return true;
   874     }
   976     }
   875 
   977 
   876     /// Finds the join node.
   978     /// Find the join node.
   877     Node findJoinNode() {
   979     inline Node findJoinNode() {
   878       Node u = _graph.source(_in_edge);
   980       Node u = _graph.source(_in_edge);
   879       Node v = _graph.target(_in_edge);
   981       Node v = _graph.target(_in_edge);
   880       while (u != v) {
   982       while (u != v) {
   881         if (_depth[u] == _depth[v]) {
   983         if (_depth[u] == _depth[v]) {
   882           u = _parent[u];
   984           u = _parent[u];
   886         else v = _parent[v];
   988         else v = _parent[v];
   887       }
   989       }
   888       return u;
   990       return u;
   889     }
   991     }
   890 
   992 
   891     /// \brief Finds the leaving edge of the cycle. Returns \c true if
   993     /// \brief Find the leaving edge of the cycle.
   892     /// the leaving edge is not the same as the entering edge.
   994     /// \return \c true if the leaving edge is not the same as the
   893     bool findLeavingEdge() {
   995     /// entering edge.
       
   996     inline bool findLeavingEdge() {
   894       // Initializing first and second nodes according to the direction
   997       // Initializing first and second nodes according to the direction
   895       // of the cycle
   998       // of the cycle
   896       if (_state[_in_edge] == STATE_LOWER) {
   999       if (_state[_in_edge] == STATE_LOWER) {
   897         first  = _graph.source(_in_edge);
  1000         first  = _graph.source(_in_edge);
   898         second = _graph.target(_in_edge);
  1001         second = _graph.target(_in_edge);
   932         }
  1035         }
   933       }
  1036       }
   934       return result;
  1037       return result;
   935     }
  1038     }
   936 
  1039 
   937     /// Changes \c flow and \c state edge maps.
  1040     /// Change \c flow and \c state edge maps.
   938     void changeFlows(bool change) {
  1041     inline void changeFlows(bool change) {
   939       // Augmenting along the cycle
  1042       // Augmenting along the cycle
   940       if (delta > 0) {
  1043       if (delta > 0) {
   941         Capacity val = _state[_in_edge] * delta;
  1044         Capacity val = _state[_in_edge] * delta;
   942         _flow[_in_edge] += val;
  1045         _flow[_in_edge] += val;
   943         for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) {
  1046         for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) {
   955       } else {
  1058       } else {
   956         _state[_in_edge] = -_state[_in_edge];
  1059         _state[_in_edge] = -_state[_in_edge];
   957       }
  1060       }
   958     }
  1061     }
   959 
  1062 
   960     /// Updates \c thread and \c parent node maps.
  1063     /// Update \c thread and \c parent node maps.
   961     void updateThreadParent() {
  1064     inline void updateThreadParent() {
   962       Node u;
  1065       Node u;
   963       v_out = _parent[u_out];
  1066       v_out = _parent[u_out];
   964 
  1067 
   965       // Handling the case when join and v_out coincide
  1068       // Handling the case when join and v_out coincide
   966       bool par_first = false;
  1069       bool par_first = false;
  1014         for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
  1117         for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
  1015         _thread[u] = right;
  1118         _thread[u] = right;
  1016       }
  1119       }
  1017     }
  1120     }
  1018 
  1121 
  1019     /// Updates \c pred_edge and \c forward node maps.
  1122     /// Update \c pred_edge and \c forward node maps.
  1020     void updatePredEdge() {
  1123     inline void updatePredEdge() {
  1021       Node u = u_out, v;
  1124       Node u = u_out, v;
  1022       while (u != u_in) {
  1125       while (u != u_in) {
  1023         v = _parent[u];
  1126         v = _parent[u];
  1024         _pred_edge[u] = _pred_edge[v];
  1127         _pred_edge[u] = _pred_edge[v];
  1025         _forward[u] = !_forward[v];
  1128         _forward[u] = !_forward[v];
  1027       }
  1130       }
  1028       _pred_edge[u_in] = _in_edge;
  1131       _pred_edge[u_in] = _in_edge;
  1029       _forward[u_in] = (u_in == _graph.source(_in_edge));
  1132       _forward[u_in] = (u_in == _graph.source(_in_edge));
  1030     }
  1133     }
  1031 
  1134 
  1032     /// Updates \c depth and \c potential node maps.
  1135     /// Update \c depth and \c potential node maps.
  1033     void updateDepthPotential() {
  1136     inline void updateDepthPotential() {
  1034       _depth[u_in] = _depth[v_in] + 1;
  1137       _depth[u_in] = _depth[v_in] + 1;
  1035       _potential[u_in] = _forward[u_in] ?
  1138       _potential[u_in] = _forward[u_in] ?
  1036         _potential[v_in] - _cost[_pred_edge[u_in]] :
  1139         _potential[v_in] - _cost[_pred_edge[u_in]] :
  1037         _potential[v_in] + _cost[_pred_edge[u_in]];
  1140         _potential[v_in] + _cost[_pred_edge[u_in]];
  1038 
  1141 
  1047         if (_depth[u] <= _depth[v_in]) break;
  1150         if (_depth[u] <= _depth[v_in]) break;
  1048         u = _thread[u];
  1151         u = _thread[u];
  1049       }
  1152       }
  1050     }
  1153     }
  1051 
  1154 
  1052     /// Executes the algorithm.
  1155     /// Execute the algorithm.
  1053     bool start(PivotRuleEnum pivot_rule) {
  1156     bool start(PivotRuleEnum pivot_rule) {
       
  1157       // Selecting the pivot rule implementation
  1054       switch (pivot_rule) {
  1158       switch (pivot_rule) {
  1055         case FIRST_ELIGIBLE_PIVOT:
  1159         case FIRST_ELIGIBLE_PIVOT:
  1056           return start<FirstEligiblePivotRule>();
  1160           return start<FirstEligiblePivotRule>();
  1057         case BEST_ELIGIBLE_PIVOT:
  1161         case BEST_ELIGIBLE_PIVOT:
  1058           return start<BestEligiblePivotRule>();
  1162           return start<BestEligiblePivotRule>();
  1059         case BLOCK_SEARCH_PIVOT:
  1163         case BLOCK_SEARCH_PIVOT:
  1060           return start<BlockSearchPivotRule>();
  1164           return start<BlockSearchPivotRule>();
  1061         case LIMITED_SEARCH_PIVOT:
       
  1062           return start<LimitedSearchPivotRule>();
       
  1063         case CANDIDATE_LIST_PIVOT:
  1165         case CANDIDATE_LIST_PIVOT:
  1064           return start<CandidateListPivotRule>();
  1166           return start<CandidateListPivotRule>();
  1065         case COMBINED_PIVOT:
  1167         case ALTERING_LIST_PIVOT:
  1066           if ( countEdges(_graph) / countNodes(_graph) <=
  1168           return start<AlteringListPivotRule>();
  1067                COMBINED_PIVOT_MAX_DEG )
       
  1068             return start<LimitedSearchPivotRule>();
       
  1069           else
       
  1070             return start<BlockSearchPivotRule>();
       
  1071       }
  1169       }
  1072       return false;
  1170       return false;
  1073     }
  1171     }
  1074 
  1172 
  1075     template<class PivotRuleImplementation>
  1173     template<class PivotRuleImplementation>
  1076     bool start() {
  1174     bool start() {
  1077       PivotRuleImplementation pivot(*this);
  1175       PivotRuleImplementation pivot(*this, _edges);
  1078 
  1176 
  1079       // Executing the network simplex algorithm
  1177       // Executing the network simplex algorithm
  1080       while (pivot.findEnteringEdge()) {
  1178       while (pivot.findEnteringEdge()) {
  1081         join = findJoinNode();
  1179         join = findJoinNode();
  1082         bool change = findLeavingEdge();
  1180         bool change = findLeavingEdge();