lemon/network_simplex.h
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 24 Apr 2009 12:23:53 +0200
changeset 615 e3d9bff447ed
parent 609 e6927fe719e6
child 618 b95898314e09
permissions -rw-r--r--
Exploit the changes of #190 in MCF test file (#234, #190)
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2009
     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 Network Simplex algorithm for finding a minimum cost flow.
    26 
    27 #include <vector>
    28 #include <limits>
    29 #include <algorithm>
    30 
    31 #include <lemon/core.h>
    32 #include <lemon/math.h>
    33 #include <lemon/maps.h>
    34 #include <lemon/circulation.h>
    35 #include <lemon/adaptors.h>
    36 
    37 namespace lemon {
    38 
    39   /// \addtogroup min_cost_flow
    40   /// @{
    41 
    42   /// \brief Implementation of the primal Network Simplex algorithm
    43   /// for finding a \ref min_cost_flow "minimum cost flow".
    44   ///
    45   /// \ref NetworkSimplex implements the primal Network Simplex algorithm
    46   /// for finding a \ref min_cost_flow "minimum cost flow".
    47   /// This algorithm is a specialized version of the linear programming
    48   /// simplex method directly for the minimum cost flow problem.
    49   /// It is one of the most efficient solution methods.
    50   ///
    51   /// In general this class is the fastest implementation available
    52   /// in LEMON for the minimum cost flow problem.
    53   /// Moreover it supports both direction of the supply/demand inequality
    54   /// constraints. For more information see \ref ProblemType.
    55   ///
    56   /// \tparam GR The digraph type the algorithm runs on.
    57   /// \tparam F The value type used for flow amounts, capacity bounds
    58   /// and supply values in the algorithm. By default it is \c int.
    59   /// \tparam C The value type used for costs and potentials in the
    60   /// algorithm. By default it is the same as \c F.
    61   ///
    62   /// \warning Both value types must be signed and all input data must
    63   /// be integer.
    64   ///
    65   /// \note %NetworkSimplex provides five different pivot rule
    66   /// implementations, from which the most efficient one is used
    67   /// by default. For more information see \ref PivotRule.
    68   template <typename GR, typename F = int, typename C = F>
    69   class NetworkSimplex
    70   {
    71   public:
    72 
    73     /// The flow type of the algorithm
    74     typedef F Flow;
    75     /// The cost type of the algorithm
    76     typedef C Cost;
    77 #ifdef DOXYGEN
    78     /// The type of the flow map
    79     typedef GR::ArcMap<Flow> FlowMap;
    80     /// The type of the potential map
    81     typedef GR::NodeMap<Cost> PotentialMap;
    82 #else
    83     /// The type of the flow map
    84     typedef typename GR::template ArcMap<Flow> FlowMap;
    85     /// The type of the potential map
    86     typedef typename GR::template NodeMap<Cost> PotentialMap;
    87 #endif
    88 
    89   public:
    90 
    91     /// \brief Enum type for selecting the pivot rule.
    92     ///
    93     /// Enum type for selecting the pivot rule for the \ref run()
    94     /// function.
    95     ///
    96     /// \ref NetworkSimplex provides five different pivot rule
    97     /// implementations that significantly affect the running time
    98     /// of the algorithm.
    99     /// By default \ref BLOCK_SEARCH "Block Search" is used, which
   100     /// proved to be the most efficient and the most robust on various
   101     /// test inputs according to our benchmark tests.
   102     /// However another pivot rule can be selected using the \ref run()
   103     /// function with the proper parameter.
   104     enum PivotRule {
   105 
   106       /// The First Eligible pivot rule.
   107       /// The next eligible arc is selected in a wraparound fashion
   108       /// in every iteration.
   109       FIRST_ELIGIBLE,
   110 
   111       /// The Best Eligible pivot rule.
   112       /// The best eligible arc is selected in every iteration.
   113       BEST_ELIGIBLE,
   114 
   115       /// The Block Search pivot rule.
   116       /// A specified number of arcs are examined in every iteration
   117       /// in a wraparound fashion and the best eligible arc is selected
   118       /// from this block.
   119       BLOCK_SEARCH,
   120 
   121       /// The Candidate List pivot rule.
   122       /// In a major iteration a candidate list is built from eligible arcs
   123       /// in a wraparound fashion and in the following minor iterations
   124       /// the best eligible arc is selected from this list.
   125       CANDIDATE_LIST,
   126 
   127       /// The Altering Candidate List pivot rule.
   128       /// It is a modified version of the Candidate List method.
   129       /// It keeps only the several best eligible arcs from the former
   130       /// candidate list and extends this list in every iteration.
   131       ALTERING_LIST
   132     };
   133     
   134     /// \brief Enum type for selecting the problem type.
   135     ///
   136     /// Enum type for selecting the problem type, i.e. the direction of
   137     /// the inequalities in the supply/demand constraints of the
   138     /// \ref min_cost_flow "minimum cost flow problem".
   139     ///
   140     /// The default problem type is \c GEQ, since this form is supported
   141     /// by other minimum cost flow algorithms and the \ref Circulation
   142     /// algorithm as well.
   143     /// The \c LEQ problem type can be selected using the \ref problemType()
   144     /// function.
   145     ///
   146     /// Note that the equality form is a special case of both problem type.
   147     enum ProblemType {
   148 
   149       /// This option means that there are "<em>greater or equal</em>"
   150       /// constraints in the defintion, i.e. the exact formulation of the
   151       /// problem is the following.
   152       /**
   153           \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
   154           \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
   155               sup(u) \quad \forall u\in V \f]
   156           \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
   157       */
   158       /// It means that the total demand must be greater or equal to the 
   159       /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
   160       /// negative) and all the supplies have to be carried out from 
   161       /// the supply nodes, but there could be demands that are not 
   162       /// satisfied.
   163       GEQ,
   164       /// It is just an alias for the \c GEQ option.
   165       CARRY_SUPPLIES = GEQ,
   166 
   167       /// This option means that there are "<em>less or equal</em>"
   168       /// constraints in the defintion, i.e. the exact formulation of the
   169       /// problem is the following.
   170       /**
   171           \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
   172           \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
   173               sup(u) \quad \forall u\in V \f]
   174           \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
   175       */
   176       /// It means that the total demand must be less or equal to the 
   177       /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
   178       /// positive) and all the demands have to be satisfied, but there
   179       /// could be supplies that are not carried out from the supply
   180       /// nodes.
   181       LEQ,
   182       /// It is just an alias for the \c LEQ option.
   183       SATISFY_DEMANDS = LEQ
   184     };
   185 
   186   private:
   187 
   188     TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   189 
   190     typedef typename GR::template ArcMap<Flow> FlowArcMap;
   191     typedef typename GR::template ArcMap<Cost> CostArcMap;
   192     typedef typename GR::template NodeMap<Flow> FlowNodeMap;
   193 
   194     typedef std::vector<Arc> ArcVector;
   195     typedef std::vector<Node> NodeVector;
   196     typedef std::vector<int> IntVector;
   197     typedef std::vector<bool> BoolVector;
   198     typedef std::vector<Flow> FlowVector;
   199     typedef std::vector<Cost> CostVector;
   200 
   201     // State constants for arcs
   202     enum ArcStateEnum {
   203       STATE_UPPER = -1,
   204       STATE_TREE  =  0,
   205       STATE_LOWER =  1
   206     };
   207 
   208   private:
   209 
   210     // Data related to the underlying digraph
   211     const GR &_graph;
   212     int _node_num;
   213     int _arc_num;
   214 
   215     // Parameters of the problem
   216     FlowArcMap *_plower;
   217     FlowArcMap *_pupper;
   218     CostArcMap *_pcost;
   219     FlowNodeMap *_psupply;
   220     bool _pstsup;
   221     Node _psource, _ptarget;
   222     Flow _pstflow;
   223     ProblemType _ptype;
   224 
   225     // Result maps
   226     FlowMap *_flow_map;
   227     PotentialMap *_potential_map;
   228     bool _local_flow;
   229     bool _local_potential;
   230 
   231     // Data structures for storing the digraph
   232     IntNodeMap _node_id;
   233     ArcVector _arc_ref;
   234     IntVector _source;
   235     IntVector _target;
   236 
   237     // Node and arc data
   238     FlowVector _cap;
   239     CostVector _cost;
   240     FlowVector _supply;
   241     FlowVector _flow;
   242     CostVector _pi;
   243 
   244     // Data for storing the spanning tree structure
   245     IntVector _parent;
   246     IntVector _pred;
   247     IntVector _thread;
   248     IntVector _rev_thread;
   249     IntVector _succ_num;
   250     IntVector _last_succ;
   251     IntVector _dirty_revs;
   252     BoolVector _forward;
   253     IntVector _state;
   254     int _root;
   255 
   256     // Temporary data used in the current pivot iteration
   257     int in_arc, join, u_in, v_in, u_out, v_out;
   258     int first, second, right, last;
   259     int stem, par_stem, new_stem;
   260     Flow delta;
   261 
   262   private:
   263 
   264     // Implementation of the First Eligible pivot rule
   265     class FirstEligiblePivotRule
   266     {
   267     private:
   268 
   269       // References to the NetworkSimplex class
   270       const IntVector  &_source;
   271       const IntVector  &_target;
   272       const CostVector &_cost;
   273       const IntVector  &_state;
   274       const CostVector &_pi;
   275       int &_in_arc;
   276       int _arc_num;
   277 
   278       // Pivot rule data
   279       int _next_arc;
   280 
   281     public:
   282 
   283       // Constructor
   284       FirstEligiblePivotRule(NetworkSimplex &ns) :
   285         _source(ns._source), _target(ns._target),
   286         _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   287         _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
   288       {}
   289 
   290       // Find next entering arc
   291       bool findEnteringArc() {
   292         Cost c;
   293         for (int e = _next_arc; e < _arc_num; ++e) {
   294           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   295           if (c < 0) {
   296             _in_arc = e;
   297             _next_arc = e + 1;
   298             return true;
   299           }
   300         }
   301         for (int e = 0; e < _next_arc; ++e) {
   302           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   303           if (c < 0) {
   304             _in_arc = e;
   305             _next_arc = e + 1;
   306             return true;
   307           }
   308         }
   309         return false;
   310       }
   311 
   312     }; //class FirstEligiblePivotRule
   313 
   314 
   315     // Implementation of the Best Eligible pivot rule
   316     class BestEligiblePivotRule
   317     {
   318     private:
   319 
   320       // References to the NetworkSimplex class
   321       const IntVector  &_source;
   322       const IntVector  &_target;
   323       const CostVector &_cost;
   324       const IntVector  &_state;
   325       const CostVector &_pi;
   326       int &_in_arc;
   327       int _arc_num;
   328 
   329     public:
   330 
   331       // Constructor
   332       BestEligiblePivotRule(NetworkSimplex &ns) :
   333         _source(ns._source), _target(ns._target),
   334         _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   335         _in_arc(ns.in_arc), _arc_num(ns._arc_num)
   336       {}
   337 
   338       // Find next entering arc
   339       bool findEnteringArc() {
   340         Cost c, min = 0;
   341         for (int e = 0; e < _arc_num; ++e) {
   342           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   343           if (c < min) {
   344             min = c;
   345             _in_arc = e;
   346           }
   347         }
   348         return min < 0;
   349       }
   350 
   351     }; //class BestEligiblePivotRule
   352 
   353 
   354     // Implementation of the Block Search pivot rule
   355     class BlockSearchPivotRule
   356     {
   357     private:
   358 
   359       // References to the NetworkSimplex class
   360       const IntVector  &_source;
   361       const IntVector  &_target;
   362       const CostVector &_cost;
   363       const IntVector  &_state;
   364       const CostVector &_pi;
   365       int &_in_arc;
   366       int _arc_num;
   367 
   368       // Pivot rule data
   369       int _block_size;
   370       int _next_arc;
   371 
   372     public:
   373 
   374       // Constructor
   375       BlockSearchPivotRule(NetworkSimplex &ns) :
   376         _source(ns._source), _target(ns._target),
   377         _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   378         _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
   379       {
   380         // The main parameters of the pivot rule
   381         const double BLOCK_SIZE_FACTOR = 2.0;
   382         const int MIN_BLOCK_SIZE = 10;
   383 
   384         _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
   385                                 MIN_BLOCK_SIZE );
   386       }
   387 
   388       // Find next entering arc
   389       bool findEnteringArc() {
   390         Cost c, min = 0;
   391         int cnt = _block_size;
   392         int e, min_arc = _next_arc;
   393         for (e = _next_arc; e < _arc_num; ++e) {
   394           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   395           if (c < min) {
   396             min = c;
   397             min_arc = e;
   398           }
   399           if (--cnt == 0) {
   400             if (min < 0) break;
   401             cnt = _block_size;
   402           }
   403         }
   404         if (min == 0 || cnt > 0) {
   405           for (e = 0; e < _next_arc; ++e) {
   406             c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   407             if (c < min) {
   408               min = c;
   409               min_arc = e;
   410             }
   411             if (--cnt == 0) {
   412               if (min < 0) break;
   413               cnt = _block_size;
   414             }
   415           }
   416         }
   417         if (min >= 0) return false;
   418         _in_arc = min_arc;
   419         _next_arc = e;
   420         return true;
   421       }
   422 
   423     }; //class BlockSearchPivotRule
   424 
   425 
   426     // Implementation of the Candidate List pivot rule
   427     class CandidateListPivotRule
   428     {
   429     private:
   430 
   431       // References to the NetworkSimplex class
   432       const IntVector  &_source;
   433       const IntVector  &_target;
   434       const CostVector &_cost;
   435       const IntVector  &_state;
   436       const CostVector &_pi;
   437       int &_in_arc;
   438       int _arc_num;
   439 
   440       // Pivot rule data
   441       IntVector _candidates;
   442       int _list_length, _minor_limit;
   443       int _curr_length, _minor_count;
   444       int _next_arc;
   445 
   446     public:
   447 
   448       /// Constructor
   449       CandidateListPivotRule(NetworkSimplex &ns) :
   450         _source(ns._source), _target(ns._target),
   451         _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   452         _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
   453       {
   454         // The main parameters of the pivot rule
   455         const double LIST_LENGTH_FACTOR = 1.0;
   456         const int MIN_LIST_LENGTH = 10;
   457         const double MINOR_LIMIT_FACTOR = 0.1;
   458         const int MIN_MINOR_LIMIT = 3;
   459 
   460         _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_arc_num)),
   461                                  MIN_LIST_LENGTH );
   462         _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
   463                                  MIN_MINOR_LIMIT );
   464         _curr_length = _minor_count = 0;
   465         _candidates.resize(_list_length);
   466       }
   467 
   468       /// Find next entering arc
   469       bool findEnteringArc() {
   470         Cost min, c;
   471         int e, min_arc = _next_arc;
   472         if (_curr_length > 0 && _minor_count < _minor_limit) {
   473           // Minor iteration: select the best eligible arc from the
   474           // current candidate list
   475           ++_minor_count;
   476           min = 0;
   477           for (int i = 0; i < _curr_length; ++i) {
   478             e = _candidates[i];
   479             c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   480             if (c < min) {
   481               min = c;
   482               min_arc = e;
   483             }
   484             if (c >= 0) {
   485               _candidates[i--] = _candidates[--_curr_length];
   486             }
   487           }
   488           if (min < 0) {
   489             _in_arc = min_arc;
   490             return true;
   491           }
   492         }
   493 
   494         // Major iteration: build a new candidate list
   495         min = 0;
   496         _curr_length = 0;
   497         for (e = _next_arc; e < _arc_num; ++e) {
   498           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   499           if (c < 0) {
   500             _candidates[_curr_length++] = e;
   501             if (c < min) {
   502               min = c;
   503               min_arc = e;
   504             }
   505             if (_curr_length == _list_length) break;
   506           }
   507         }
   508         if (_curr_length < _list_length) {
   509           for (e = 0; e < _next_arc; ++e) {
   510             c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   511             if (c < 0) {
   512               _candidates[_curr_length++] = e;
   513               if (c < min) {
   514                 min = c;
   515                 min_arc = e;
   516               }
   517               if (_curr_length == _list_length) break;
   518             }
   519           }
   520         }
   521         if (_curr_length == 0) return false;
   522         _minor_count = 1;
   523         _in_arc = min_arc;
   524         _next_arc = e;
   525         return true;
   526       }
   527 
   528     }; //class CandidateListPivotRule
   529 
   530 
   531     // Implementation of the Altering Candidate List pivot rule
   532     class AlteringListPivotRule
   533     {
   534     private:
   535 
   536       // References to the NetworkSimplex class
   537       const IntVector  &_source;
   538       const IntVector  &_target;
   539       const CostVector &_cost;
   540       const IntVector  &_state;
   541       const CostVector &_pi;
   542       int &_in_arc;
   543       int _arc_num;
   544 
   545       // Pivot rule data
   546       int _block_size, _head_length, _curr_length;
   547       int _next_arc;
   548       IntVector _candidates;
   549       CostVector _cand_cost;
   550 
   551       // Functor class to compare arcs during sort of the candidate list
   552       class SortFunc
   553       {
   554       private:
   555         const CostVector &_map;
   556       public:
   557         SortFunc(const CostVector &map) : _map(map) {}
   558         bool operator()(int left, int right) {
   559           return _map[left] > _map[right];
   560         }
   561       };
   562 
   563       SortFunc _sort_func;
   564 
   565     public:
   566 
   567       // Constructor
   568       AlteringListPivotRule(NetworkSimplex &ns) :
   569         _source(ns._source), _target(ns._target),
   570         _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   571         _in_arc(ns.in_arc), _arc_num(ns._arc_num),
   572         _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
   573       {
   574         // The main parameters of the pivot rule
   575         const double BLOCK_SIZE_FACTOR = 1.5;
   576         const int MIN_BLOCK_SIZE = 10;
   577         const double HEAD_LENGTH_FACTOR = 0.1;
   578         const int MIN_HEAD_LENGTH = 3;
   579 
   580         _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
   581                                 MIN_BLOCK_SIZE );
   582         _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
   583                                  MIN_HEAD_LENGTH );
   584         _candidates.resize(_head_length + _block_size);
   585         _curr_length = 0;
   586       }
   587 
   588       // Find next entering arc
   589       bool findEnteringArc() {
   590         // Check the current candidate list
   591         int e;
   592         for (int i = 0; i < _curr_length; ++i) {
   593           e = _candidates[i];
   594           _cand_cost[e] = _state[e] *
   595             (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   596           if (_cand_cost[e] >= 0) {
   597             _candidates[i--] = _candidates[--_curr_length];
   598           }
   599         }
   600 
   601         // Extend the list
   602         int cnt = _block_size;
   603         int last_arc = 0;
   604         int limit = _head_length;
   605 
   606         for (int e = _next_arc; e < _arc_num; ++e) {
   607           _cand_cost[e] = _state[e] *
   608             (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   609           if (_cand_cost[e] < 0) {
   610             _candidates[_curr_length++] = e;
   611             last_arc = e;
   612           }
   613           if (--cnt == 0) {
   614             if (_curr_length > limit) break;
   615             limit = 0;
   616             cnt = _block_size;
   617           }
   618         }
   619         if (_curr_length <= limit) {
   620           for (int e = 0; e < _next_arc; ++e) {
   621             _cand_cost[e] = _state[e] *
   622               (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   623             if (_cand_cost[e] < 0) {
   624               _candidates[_curr_length++] = e;
   625               last_arc = e;
   626             }
   627             if (--cnt == 0) {
   628               if (_curr_length > limit) break;
   629               limit = 0;
   630               cnt = _block_size;
   631             }
   632           }
   633         }
   634         if (_curr_length == 0) return false;
   635         _next_arc = last_arc + 1;
   636 
   637         // Make heap of the candidate list (approximating a partial sort)
   638         make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
   639                    _sort_func );
   640 
   641         // Pop the first element of the heap
   642         _in_arc = _candidates[0];
   643         pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
   644                   _sort_func );
   645         _curr_length = std::min(_head_length, _curr_length - 1);
   646         return true;
   647       }
   648 
   649     }; //class AlteringListPivotRule
   650 
   651   public:
   652 
   653     /// \brief Constructor.
   654     ///
   655     /// The constructor of the class.
   656     ///
   657     /// \param graph The digraph the algorithm runs on.
   658     NetworkSimplex(const GR& graph) :
   659       _graph(graph),
   660       _plower(NULL), _pupper(NULL), _pcost(NULL),
   661       _psupply(NULL), _pstsup(false), _ptype(GEQ),
   662       _flow_map(NULL), _potential_map(NULL),
   663       _local_flow(false), _local_potential(false),
   664       _node_id(graph)
   665     {
   666       LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
   667                    std::numeric_limits<Flow>::is_signed,
   668         "The flow type of NetworkSimplex must be signed integer");
   669       LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
   670                    std::numeric_limits<Cost>::is_signed,
   671         "The cost type of NetworkSimplex must be signed integer");
   672     }
   673 
   674     /// Destructor.
   675     ~NetworkSimplex() {
   676       if (_local_flow) delete _flow_map;
   677       if (_local_potential) delete _potential_map;
   678     }
   679 
   680     /// \name Parameters
   681     /// The parameters of the algorithm can be specified using these
   682     /// functions.
   683 
   684     /// @{
   685 
   686     /// \brief Set the lower bounds on the arcs.
   687     ///
   688     /// This function sets the lower bounds on the arcs.
   689     /// If neither this function nor \ref boundMaps() is used before
   690     /// calling \ref run(), the lower bounds will be set to zero
   691     /// on all arcs.
   692     ///
   693     /// \param map An arc map storing the lower bounds.
   694     /// Its \c Value type must be convertible to the \c Flow type
   695     /// of the algorithm.
   696     ///
   697     /// \return <tt>(*this)</tt>
   698     template <typename LOWER>
   699     NetworkSimplex& lowerMap(const LOWER& map) {
   700       delete _plower;
   701       _plower = new FlowArcMap(_graph);
   702       for (ArcIt a(_graph); a != INVALID; ++a) {
   703         (*_plower)[a] = map[a];
   704       }
   705       return *this;
   706     }
   707 
   708     /// \brief Set the upper bounds (capacities) on the arcs.
   709     ///
   710     /// This function sets the upper bounds (capacities) on the arcs.
   711     /// If none of the functions \ref upperMap(), \ref capacityMap()
   712     /// and \ref boundMaps() is used before calling \ref run(),
   713     /// the upper bounds (capacities) will be set to
   714     /// \c std::numeric_limits<Flow>::max() on all arcs.
   715     ///
   716     /// \param map An arc map storing the upper bounds.
   717     /// Its \c Value type must be convertible to the \c Flow type
   718     /// of the algorithm.
   719     ///
   720     /// \return <tt>(*this)</tt>
   721     template<typename UPPER>
   722     NetworkSimplex& upperMap(const UPPER& map) {
   723       delete _pupper;
   724       _pupper = new FlowArcMap(_graph);
   725       for (ArcIt a(_graph); a != INVALID; ++a) {
   726         (*_pupper)[a] = map[a];
   727       }
   728       return *this;
   729     }
   730 
   731     /// \brief Set the upper bounds (capacities) on the arcs.
   732     ///
   733     /// This function sets the upper bounds (capacities) on the arcs.
   734     /// It is just an alias for \ref upperMap().
   735     ///
   736     /// \return <tt>(*this)</tt>
   737     template<typename CAP>
   738     NetworkSimplex& capacityMap(const CAP& map) {
   739       return upperMap(map);
   740     }
   741 
   742     /// \brief Set the lower and upper bounds on the arcs.
   743     ///
   744     /// This function sets the lower and upper bounds on the arcs.
   745     /// If neither this function nor \ref lowerMap() is used before
   746     /// calling \ref run(), the lower bounds will be set to zero
   747     /// on all arcs.
   748     /// If none of the functions \ref upperMap(), \ref capacityMap()
   749     /// and \ref boundMaps() is used before calling \ref run(),
   750     /// the upper bounds (capacities) will be set to
   751     /// \c std::numeric_limits<Flow>::max() on all arcs.
   752     ///
   753     /// \param lower An arc map storing the lower bounds.
   754     /// \param upper An arc map storing the upper bounds.
   755     ///
   756     /// The \c Value type of the maps must be convertible to the
   757     /// \c Flow type of the algorithm.
   758     ///
   759     /// \note This function is just a shortcut of calling \ref lowerMap()
   760     /// and \ref upperMap() separately.
   761     ///
   762     /// \return <tt>(*this)</tt>
   763     template <typename LOWER, typename UPPER>
   764     NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
   765       return lowerMap(lower).upperMap(upper);
   766     }
   767 
   768     /// \brief Set the costs of the arcs.
   769     ///
   770     /// This function sets the costs of the arcs.
   771     /// If it is not used before calling \ref run(), the costs
   772     /// will be set to \c 1 on all arcs.
   773     ///
   774     /// \param map An arc map storing the costs.
   775     /// Its \c Value type must be convertible to the \c Cost type
   776     /// of the algorithm.
   777     ///
   778     /// \return <tt>(*this)</tt>
   779     template<typename COST>
   780     NetworkSimplex& costMap(const COST& map) {
   781       delete _pcost;
   782       _pcost = new CostArcMap(_graph);
   783       for (ArcIt a(_graph); a != INVALID; ++a) {
   784         (*_pcost)[a] = map[a];
   785       }
   786       return *this;
   787     }
   788 
   789     /// \brief Set the supply values of the nodes.
   790     ///
   791     /// This function sets the supply values of the nodes.
   792     /// If neither this function nor \ref stSupply() is used before
   793     /// calling \ref run(), the supply of each node will be set to zero.
   794     /// (It makes sense only if non-zero lower bounds are given.)
   795     ///
   796     /// \param map A node map storing the supply values.
   797     /// Its \c Value type must be convertible to the \c Flow type
   798     /// of the algorithm.
   799     ///
   800     /// \return <tt>(*this)</tt>
   801     template<typename SUP>
   802     NetworkSimplex& supplyMap(const SUP& map) {
   803       delete _psupply;
   804       _pstsup = false;
   805       _psupply = new FlowNodeMap(_graph);
   806       for (NodeIt n(_graph); n != INVALID; ++n) {
   807         (*_psupply)[n] = map[n];
   808       }
   809       return *this;
   810     }
   811 
   812     /// \brief Set single source and target nodes and a supply value.
   813     ///
   814     /// This function sets a single source node and a single target node
   815     /// and the required flow value.
   816     /// If neither this function nor \ref supplyMap() is used before
   817     /// calling \ref run(), the supply of each node will be set to zero.
   818     /// (It makes sense only if non-zero lower bounds are given.)
   819     ///
   820     /// \param s The source node.
   821     /// \param t The target node.
   822     /// \param k The required amount of flow from node \c s to node \c t
   823     /// (i.e. the supply of \c s and the demand of \c t).
   824     ///
   825     /// \return <tt>(*this)</tt>
   826     NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) {
   827       delete _psupply;
   828       _psupply = NULL;
   829       _pstsup = true;
   830       _psource = s;
   831       _ptarget = t;
   832       _pstflow = k;
   833       return *this;
   834     }
   835     
   836     /// \brief Set the problem type.
   837     ///
   838     /// This function sets the problem type for the algorithm.
   839     /// If it is not used before calling \ref run(), the \ref GEQ problem
   840     /// type will be used.
   841     ///
   842     /// For more information see \ref ProblemType.
   843     ///
   844     /// \return <tt>(*this)</tt>
   845     NetworkSimplex& problemType(ProblemType problem_type) {
   846       _ptype = problem_type;
   847       return *this;
   848     }
   849 
   850     /// \brief Set the flow map.
   851     ///
   852     /// This function sets the flow map.
   853     /// If it is not used before calling \ref run(), an instance will
   854     /// be allocated automatically. The destructor deallocates this
   855     /// automatically allocated map, of course.
   856     ///
   857     /// \return <tt>(*this)</tt>
   858     NetworkSimplex& flowMap(FlowMap& map) {
   859       if (_local_flow) {
   860         delete _flow_map;
   861         _local_flow = false;
   862       }
   863       _flow_map = &map;
   864       return *this;
   865     }
   866 
   867     /// \brief Set the potential map.
   868     ///
   869     /// This function sets the potential map, which is used for storing
   870     /// the dual solution.
   871     /// If it is not used before calling \ref run(), an instance will
   872     /// be allocated automatically. The destructor deallocates this
   873     /// automatically allocated map, of course.
   874     ///
   875     /// \return <tt>(*this)</tt>
   876     NetworkSimplex& potentialMap(PotentialMap& map) {
   877       if (_local_potential) {
   878         delete _potential_map;
   879         _local_potential = false;
   880       }
   881       _potential_map = &map;
   882       return *this;
   883     }
   884     
   885     /// @}
   886 
   887     /// \name Execution Control
   888     /// The algorithm can be executed using \ref run().
   889 
   890     /// @{
   891 
   892     /// \brief Run the algorithm.
   893     ///
   894     /// This function runs the algorithm.
   895     /// The paramters can be specified using functions \ref lowerMap(),
   896     /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
   897     /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 
   898     /// \ref problemType(), \ref flowMap() and \ref potentialMap().
   899     /// For example,
   900     /// \code
   901     ///   NetworkSimplex<ListDigraph> ns(graph);
   902     ///   ns.boundMaps(lower, upper).costMap(cost)
   903     ///     .supplyMap(sup).run();
   904     /// \endcode
   905     ///
   906     /// This function can be called more than once. All the parameters
   907     /// that have been given are kept for the next call, unless
   908     /// \ref reset() is called, thus only the modified parameters
   909     /// have to be set again. See \ref reset() for examples.
   910     ///
   911     /// \param pivot_rule The pivot rule that will be used during the
   912     /// algorithm. For more information see \ref PivotRule.
   913     ///
   914     /// \return \c true if a feasible flow can be found.
   915     bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
   916       return init() && start(pivot_rule);
   917     }
   918 
   919     /// \brief Reset all the parameters that have been given before.
   920     ///
   921     /// This function resets all the paramaters that have been given
   922     /// before using functions \ref lowerMap(), \ref upperMap(),
   923     /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
   924     /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 
   925     /// \ref flowMap() and \ref potentialMap().
   926     ///
   927     /// It is useful for multiple run() calls. If this function is not
   928     /// used, all the parameters given before are kept for the next
   929     /// \ref run() call.
   930     ///
   931     /// For example,
   932     /// \code
   933     ///   NetworkSimplex<ListDigraph> ns(graph);
   934     ///
   935     ///   // First run
   936     ///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)
   937     ///     .supplyMap(sup).run();
   938     ///
   939     ///   // Run again with modified cost map (reset() is not called,
   940     ///   // so only the cost map have to be set again)
   941     ///   cost[e] += 100;
   942     ///   ns.costMap(cost).run();
   943     ///
   944     ///   // Run again from scratch using reset()
   945     ///   // (the lower bounds will be set to zero on all arcs)
   946     ///   ns.reset();
   947     ///   ns.capacityMap(cap).costMap(cost)
   948     ///     .supplyMap(sup).run();
   949     /// \endcode
   950     ///
   951     /// \return <tt>(*this)</tt>
   952     NetworkSimplex& reset() {
   953       delete _plower;
   954       delete _pupper;
   955       delete _pcost;
   956       delete _psupply;
   957       _plower = NULL;
   958       _pupper = NULL;
   959       _pcost = NULL;
   960       _psupply = NULL;
   961       _pstsup = false;
   962       _ptype = GEQ;
   963       if (_local_flow) delete _flow_map;
   964       if (_local_potential) delete _potential_map;
   965       _flow_map = NULL;
   966       _potential_map = NULL;
   967       _local_flow = false;
   968       _local_potential = false;
   969 
   970       return *this;
   971     }
   972 
   973     /// @}
   974 
   975     /// \name Query Functions
   976     /// The results of the algorithm can be obtained using these
   977     /// functions.\n
   978     /// The \ref run() function must be called before using them.
   979 
   980     /// @{
   981 
   982     /// \brief Return the total cost of the found flow.
   983     ///
   984     /// This function returns the total cost of the found flow.
   985     /// The complexity of the function is O(e).
   986     ///
   987     /// \note The return type of the function can be specified as a
   988     /// template parameter. For example,
   989     /// \code
   990     ///   ns.totalCost<double>();
   991     /// \endcode
   992     /// It is useful if the total cost cannot be stored in the \c Cost
   993     /// type of the algorithm, which is the default return type of the
   994     /// function.
   995     ///
   996     /// \pre \ref run() must be called before using this function.
   997     template <typename Num>
   998     Num totalCost() const {
   999       Num c = 0;
  1000       if (_pcost) {
  1001         for (ArcIt e(_graph); e != INVALID; ++e)
  1002           c += (*_flow_map)[e] * (*_pcost)[e];
  1003       } else {
  1004         for (ArcIt e(_graph); e != INVALID; ++e)
  1005           c += (*_flow_map)[e];
  1006       }
  1007       return c;
  1008     }
  1009 
  1010 #ifndef DOXYGEN
  1011     Cost totalCost() const {
  1012       return totalCost<Cost>();
  1013     }
  1014 #endif
  1015 
  1016     /// \brief Return the flow on the given arc.
  1017     ///
  1018     /// This function returns the flow on the given arc.
  1019     ///
  1020     /// \pre \ref run() must be called before using this function.
  1021     Flow flow(const Arc& a) const {
  1022       return (*_flow_map)[a];
  1023     }
  1024 
  1025     /// \brief Return a const reference to the flow map.
  1026     ///
  1027     /// This function returns a const reference to an arc map storing
  1028     /// the found flow.
  1029     ///
  1030     /// \pre \ref run() must be called before using this function.
  1031     const FlowMap& flowMap() const {
  1032       return *_flow_map;
  1033     }
  1034 
  1035     /// \brief Return the potential (dual value) of the given node.
  1036     ///
  1037     /// This function returns the potential (dual value) of the
  1038     /// given node.
  1039     ///
  1040     /// \pre \ref run() must be called before using this function.
  1041     Cost potential(const Node& n) const {
  1042       return (*_potential_map)[n];
  1043     }
  1044 
  1045     /// \brief Return a const reference to the potential map
  1046     /// (the dual solution).
  1047     ///
  1048     /// This function returns a const reference to a node map storing
  1049     /// the found potentials, which form the dual solution of the
  1050     /// \ref min_cost_flow "minimum cost flow" problem.
  1051     ///
  1052     /// \pre \ref run() must be called before using this function.
  1053     const PotentialMap& potentialMap() const {
  1054       return *_potential_map;
  1055     }
  1056 
  1057     /// @}
  1058 
  1059   private:
  1060 
  1061     // Initialize internal data structures
  1062     bool init() {
  1063       // Initialize result maps
  1064       if (!_flow_map) {
  1065         _flow_map = new FlowMap(_graph);
  1066         _local_flow = true;
  1067       }
  1068       if (!_potential_map) {
  1069         _potential_map = new PotentialMap(_graph);
  1070         _local_potential = true;
  1071       }
  1072 
  1073       // Initialize vectors
  1074       _node_num = countNodes(_graph);
  1075       _arc_num = countArcs(_graph);
  1076       int all_node_num = _node_num + 1;
  1077       int all_arc_num = _arc_num + _node_num;
  1078       if (_node_num == 0) return false;
  1079 
  1080       _arc_ref.resize(_arc_num);
  1081       _source.resize(all_arc_num);
  1082       _target.resize(all_arc_num);
  1083 
  1084       _cap.resize(all_arc_num);
  1085       _cost.resize(all_arc_num);
  1086       _supply.resize(all_node_num);
  1087       _flow.resize(all_arc_num);
  1088       _pi.resize(all_node_num);
  1089 
  1090       _parent.resize(all_node_num);
  1091       _pred.resize(all_node_num);
  1092       _forward.resize(all_node_num);
  1093       _thread.resize(all_node_num);
  1094       _rev_thread.resize(all_node_num);
  1095       _succ_num.resize(all_node_num);
  1096       _last_succ.resize(all_node_num);
  1097       _state.resize(all_arc_num);
  1098 
  1099       // Initialize node related data
  1100       bool valid_supply = true;
  1101       Flow sum_supply = 0;
  1102       if (!_pstsup && !_psupply) {
  1103         _pstsup = true;
  1104         _psource = _ptarget = NodeIt(_graph);
  1105         _pstflow = 0;
  1106       }
  1107       if (_psupply) {
  1108         int i = 0;
  1109         for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
  1110           _node_id[n] = i;
  1111           _supply[i] = (*_psupply)[n];
  1112           sum_supply += _supply[i];
  1113         }
  1114         valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
  1115                        (_ptype == LEQ && sum_supply >= 0);
  1116       } else {
  1117         int i = 0;
  1118         for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
  1119           _node_id[n] = i;
  1120           _supply[i] = 0;
  1121         }
  1122         _supply[_node_id[_psource]] =  _pstflow;
  1123         _supply[_node_id[_ptarget]] = -_pstflow;
  1124       }
  1125       if (!valid_supply) return false;
  1126 
  1127       // Infinite capacity value
  1128       Flow inf_cap =
  1129         std::numeric_limits<Flow>::has_infinity ?
  1130         std::numeric_limits<Flow>::infinity() :
  1131         std::numeric_limits<Flow>::max();
  1132 
  1133       // Initialize artifical cost
  1134       Cost art_cost;
  1135       if (std::numeric_limits<Cost>::is_exact) {
  1136         art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
  1137       } else {
  1138         art_cost = std::numeric_limits<Cost>::min();
  1139         for (int i = 0; i != _arc_num; ++i) {
  1140           if (_cost[i] > art_cost) art_cost = _cost[i];
  1141         }
  1142         art_cost = (art_cost + 1) * _node_num;
  1143       }
  1144 
  1145       // Run Circulation to check if a feasible solution exists
  1146       typedef ConstMap<Arc, Flow> ConstArcMap;
  1147       ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap);
  1148       FlowNodeMap *csup = NULL;
  1149       bool local_csup = false;
  1150       if (_psupply) {
  1151         csup = _psupply;
  1152       } else {
  1153         csup = new FlowNodeMap(_graph, 0);
  1154         (*csup)[_psource] =  _pstflow;
  1155         (*csup)[_ptarget] = -_pstflow;
  1156         local_csup = true;
  1157       }
  1158       bool circ_result = false;
  1159       if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
  1160         // GEQ problem type
  1161         if (_plower) {
  1162           if (_pupper) {
  1163             Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
  1164               circ(_graph, *_plower, *_pupper, *csup);
  1165             circ_result = circ.run();
  1166           } else {
  1167             Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
  1168               circ(_graph, *_plower, inf_arc_map, *csup);
  1169             circ_result = circ.run();
  1170           }
  1171         } else {
  1172           if (_pupper) {
  1173             Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
  1174               circ(_graph, zero_arc_map, *_pupper, *csup);
  1175             circ_result = circ.run();
  1176           } else {
  1177             Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
  1178               circ(_graph, zero_arc_map, inf_arc_map, *csup);
  1179             circ_result = circ.run();
  1180           }
  1181         }
  1182       } else {
  1183         // LEQ problem type
  1184         typedef ReverseDigraph<const GR> RevGraph;
  1185         typedef NegMap<FlowNodeMap> NegNodeMap;
  1186         RevGraph rgraph(_graph);
  1187         NegNodeMap neg_csup(*csup);
  1188         if (_plower) {
  1189           if (_pupper) {
  1190             Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
  1191               circ(rgraph, *_plower, *_pupper, neg_csup);
  1192             circ_result = circ.run();
  1193           } else {
  1194             Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
  1195               circ(rgraph, *_plower, inf_arc_map, neg_csup);
  1196             circ_result = circ.run();
  1197           }
  1198         } else {
  1199           if (_pupper) {
  1200             Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
  1201               circ(rgraph, zero_arc_map, *_pupper, neg_csup);
  1202             circ_result = circ.run();
  1203           } else {
  1204             Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
  1205               circ(rgraph, zero_arc_map, inf_arc_map, neg_csup);
  1206             circ_result = circ.run();
  1207           }
  1208         }
  1209       }
  1210       if (local_csup) delete csup;
  1211       if (!circ_result) return false;
  1212 
  1213       // Set data for the artificial root node
  1214       _root = _node_num;
  1215       _parent[_root] = -1;
  1216       _pred[_root] = -1;
  1217       _thread[_root] = 0;
  1218       _rev_thread[0] = _root;
  1219       _succ_num[_root] = all_node_num;
  1220       _last_succ[_root] = _root - 1;
  1221       _supply[_root] = -sum_supply;
  1222       if (sum_supply < 0) {
  1223         _pi[_root] = -art_cost;
  1224       } else {
  1225         _pi[_root] = art_cost;
  1226       }
  1227 
  1228       // Store the arcs in a mixed order
  1229       int k = std::max(int(sqrt(_arc_num)), 10);
  1230       int i = 0;
  1231       for (ArcIt e(_graph); e != INVALID; ++e) {
  1232         _arc_ref[i] = e;
  1233         if ((i += k) >= _arc_num) i = (i % k) + 1;
  1234       }
  1235 
  1236       // Initialize arc maps
  1237       if (_pupper && _pcost) {
  1238         for (int i = 0; i != _arc_num; ++i) {
  1239           Arc e = _arc_ref[i];
  1240           _source[i] = _node_id[_graph.source(e)];
  1241           _target[i] = _node_id[_graph.target(e)];
  1242           _cap[i] = (*_pupper)[e];
  1243           _cost[i] = (*_pcost)[e];
  1244           _flow[i] = 0;
  1245           _state[i] = STATE_LOWER;
  1246         }
  1247       } else {
  1248         for (int i = 0; i != _arc_num; ++i) {
  1249           Arc e = _arc_ref[i];
  1250           _source[i] = _node_id[_graph.source(e)];
  1251           _target[i] = _node_id[_graph.target(e)];
  1252           _flow[i] = 0;
  1253           _state[i] = STATE_LOWER;
  1254         }
  1255         if (_pupper) {
  1256           for (int i = 0; i != _arc_num; ++i)
  1257             _cap[i] = (*_pupper)[_arc_ref[i]];
  1258         } else {
  1259           for (int i = 0; i != _arc_num; ++i)
  1260             _cap[i] = inf_cap;
  1261         }
  1262         if (_pcost) {
  1263           for (int i = 0; i != _arc_num; ++i)
  1264             _cost[i] = (*_pcost)[_arc_ref[i]];
  1265         } else {
  1266           for (int i = 0; i != _arc_num; ++i)
  1267             _cost[i] = 1;
  1268         }
  1269       }
  1270       
  1271       // Remove non-zero lower bounds
  1272       if (_plower) {
  1273         for (int i = 0; i != _arc_num; ++i) {
  1274           Flow c = (*_plower)[_arc_ref[i]];
  1275           if (c != 0) {
  1276             _cap[i] -= c;
  1277             _supply[_source[i]] -= c;
  1278             _supply[_target[i]] += c;
  1279           }
  1280         }
  1281       }
  1282 
  1283       // Add artificial arcs and initialize the spanning tree data structure
  1284       for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
  1285         _thread[u] = u + 1;
  1286         _rev_thread[u + 1] = u;
  1287         _succ_num[u] = 1;
  1288         _last_succ[u] = u;
  1289         _parent[u] = _root;
  1290         _pred[u] = e;
  1291         _cost[e] = art_cost;
  1292         _cap[e] = inf_cap;
  1293         _state[e] = STATE_TREE;
  1294         if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
  1295           _flow[e] = _supply[u];
  1296           _forward[u] = true;
  1297           _pi[u] = -art_cost + _pi[_root];
  1298         } else {
  1299           _flow[e] = -_supply[u];
  1300           _forward[u] = false;
  1301           _pi[u] = art_cost + _pi[_root];
  1302         }
  1303       }
  1304 
  1305       return true;
  1306     }
  1307 
  1308     // Find the join node
  1309     void findJoinNode() {
  1310       int u = _source[in_arc];
  1311       int v = _target[in_arc];
  1312       while (u != v) {
  1313         if (_succ_num[u] < _succ_num[v]) {
  1314           u = _parent[u];
  1315         } else {
  1316           v = _parent[v];
  1317         }
  1318       }
  1319       join = u;
  1320     }
  1321 
  1322     // Find the leaving arc of the cycle and returns true if the
  1323     // leaving arc is not the same as the entering arc
  1324     bool findLeavingArc() {
  1325       // Initialize first and second nodes according to the direction
  1326       // of the cycle
  1327       if (_state[in_arc] == STATE_LOWER) {
  1328         first  = _source[in_arc];
  1329         second = _target[in_arc];
  1330       } else {
  1331         first  = _target[in_arc];
  1332         second = _source[in_arc];
  1333       }
  1334       delta = _cap[in_arc];
  1335       int result = 0;
  1336       Flow d;
  1337       int e;
  1338 
  1339       // Search the cycle along the path form the first node to the root
  1340       for (int u = first; u != join; u = _parent[u]) {
  1341         e = _pred[u];
  1342         d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
  1343         if (d < delta) {
  1344           delta = d;
  1345           u_out = u;
  1346           result = 1;
  1347         }
  1348       }
  1349       // Search the cycle along the path form the second node to the root
  1350       for (int u = second; u != join; u = _parent[u]) {
  1351         e = _pred[u];
  1352         d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
  1353         if (d <= delta) {
  1354           delta = d;
  1355           u_out = u;
  1356           result = 2;
  1357         }
  1358       }
  1359 
  1360       if (result == 1) {
  1361         u_in = first;
  1362         v_in = second;
  1363       } else {
  1364         u_in = second;
  1365         v_in = first;
  1366       }
  1367       return result != 0;
  1368     }
  1369 
  1370     // Change _flow and _state vectors
  1371     void changeFlow(bool change) {
  1372       // Augment along the cycle
  1373       if (delta > 0) {
  1374         Flow val = _state[in_arc] * delta;
  1375         _flow[in_arc] += val;
  1376         for (int u = _source[in_arc]; u != join; u = _parent[u]) {
  1377           _flow[_pred[u]] += _forward[u] ? -val : val;
  1378         }
  1379         for (int u = _target[in_arc]; u != join; u = _parent[u]) {
  1380           _flow[_pred[u]] += _forward[u] ? val : -val;
  1381         }
  1382       }
  1383       // Update the state of the entering and leaving arcs
  1384       if (change) {
  1385         _state[in_arc] = STATE_TREE;
  1386         _state[_pred[u_out]] =
  1387           (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
  1388       } else {
  1389         _state[in_arc] = -_state[in_arc];
  1390       }
  1391     }
  1392 
  1393     // Update the tree structure
  1394     void updateTreeStructure() {
  1395       int u, w;
  1396       int old_rev_thread = _rev_thread[u_out];
  1397       int old_succ_num = _succ_num[u_out];
  1398       int old_last_succ = _last_succ[u_out];
  1399       v_out = _parent[u_out];
  1400 
  1401       u = _last_succ[u_in];  // the last successor of u_in
  1402       right = _thread[u];    // the node after it
  1403 
  1404       // Handle the case when old_rev_thread equals to v_in
  1405       // (it also means that join and v_out coincide)
  1406       if (old_rev_thread == v_in) {
  1407         last = _thread[_last_succ[u_out]];
  1408       } else {
  1409         last = _thread[v_in];
  1410       }
  1411 
  1412       // Update _thread and _parent along the stem nodes (i.e. the nodes
  1413       // between u_in and u_out, whose parent have to be changed)
  1414       _thread[v_in] = stem = u_in;
  1415       _dirty_revs.clear();
  1416       _dirty_revs.push_back(v_in);
  1417       par_stem = v_in;
  1418       while (stem != u_out) {
  1419         // Insert the next stem node into the thread list
  1420         new_stem = _parent[stem];
  1421         _thread[u] = new_stem;
  1422         _dirty_revs.push_back(u);
  1423 
  1424         // Remove the subtree of stem from the thread list
  1425         w = _rev_thread[stem];
  1426         _thread[w] = right;
  1427         _rev_thread[right] = w;
  1428 
  1429         // Change the parent node and shift stem nodes
  1430         _parent[stem] = par_stem;
  1431         par_stem = stem;
  1432         stem = new_stem;
  1433 
  1434         // Update u and right
  1435         u = _last_succ[stem] == _last_succ[par_stem] ?
  1436           _rev_thread[par_stem] : _last_succ[stem];
  1437         right = _thread[u];
  1438       }
  1439       _parent[u_out] = par_stem;
  1440       _thread[u] = last;
  1441       _rev_thread[last] = u;
  1442       _last_succ[u_out] = u;
  1443 
  1444       // Remove the subtree of u_out from the thread list except for
  1445       // the case when old_rev_thread equals to v_in
  1446       // (it also means that join and v_out coincide)
  1447       if (old_rev_thread != v_in) {
  1448         _thread[old_rev_thread] = right;
  1449         _rev_thread[right] = old_rev_thread;
  1450       }
  1451 
  1452       // Update _rev_thread using the new _thread values
  1453       for (int i = 0; i < int(_dirty_revs.size()); ++i) {
  1454         u = _dirty_revs[i];
  1455         _rev_thread[_thread[u]] = u;
  1456       }
  1457 
  1458       // Update _pred, _forward, _last_succ and _succ_num for the
  1459       // stem nodes from u_out to u_in
  1460       int tmp_sc = 0, tmp_ls = _last_succ[u_out];
  1461       u = u_out;
  1462       while (u != u_in) {
  1463         w = _parent[u];
  1464         _pred[u] = _pred[w];
  1465         _forward[u] = !_forward[w];
  1466         tmp_sc += _succ_num[u] - _succ_num[w];
  1467         _succ_num[u] = tmp_sc;
  1468         _last_succ[w] = tmp_ls;
  1469         u = w;
  1470       }
  1471       _pred[u_in] = in_arc;
  1472       _forward[u_in] = (u_in == _source[in_arc]);
  1473       _succ_num[u_in] = old_succ_num;
  1474 
  1475       // Set limits for updating _last_succ form v_in and v_out
  1476       // towards the root
  1477       int up_limit_in = -1;
  1478       int up_limit_out = -1;
  1479       if (_last_succ[join] == v_in) {
  1480         up_limit_out = join;
  1481       } else {
  1482         up_limit_in = join;
  1483       }
  1484 
  1485       // Update _last_succ from v_in towards the root
  1486       for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
  1487            u = _parent[u]) {
  1488         _last_succ[u] = _last_succ[u_out];
  1489       }
  1490       // Update _last_succ from v_out towards the root
  1491       if (join != old_rev_thread && v_in != old_rev_thread) {
  1492         for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
  1493              u = _parent[u]) {
  1494           _last_succ[u] = old_rev_thread;
  1495         }
  1496       } else {
  1497         for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
  1498              u = _parent[u]) {
  1499           _last_succ[u] = _last_succ[u_out];
  1500         }
  1501       }
  1502 
  1503       // Update _succ_num from v_in to join
  1504       for (u = v_in; u != join; u = _parent[u]) {
  1505         _succ_num[u] += old_succ_num;
  1506       }
  1507       // Update _succ_num from v_out to join
  1508       for (u = v_out; u != join; u = _parent[u]) {
  1509         _succ_num[u] -= old_succ_num;
  1510       }
  1511     }
  1512 
  1513     // Update potentials
  1514     void updatePotential() {
  1515       Cost sigma = _forward[u_in] ?
  1516         _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
  1517         _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
  1518       // Update potentials in the subtree, which has been moved
  1519       int end = _thread[_last_succ[u_in]];
  1520       for (int u = u_in; u != end; u = _thread[u]) {
  1521         _pi[u] += sigma;
  1522       }
  1523     }
  1524 
  1525     // Execute the algorithm
  1526     bool start(PivotRule pivot_rule) {
  1527       // Select the pivot rule implementation
  1528       switch (pivot_rule) {
  1529         case FIRST_ELIGIBLE:
  1530           return start<FirstEligiblePivotRule>();
  1531         case BEST_ELIGIBLE:
  1532           return start<BestEligiblePivotRule>();
  1533         case BLOCK_SEARCH:
  1534           return start<BlockSearchPivotRule>();
  1535         case CANDIDATE_LIST:
  1536           return start<CandidateListPivotRule>();
  1537         case ALTERING_LIST:
  1538           return start<AlteringListPivotRule>();
  1539       }
  1540       return false;
  1541     }
  1542 
  1543     template <typename PivotRuleImpl>
  1544     bool start() {
  1545       PivotRuleImpl pivot(*this);
  1546 
  1547       // Execute the Network Simplex algorithm
  1548       while (pivot.findEnteringArc()) {
  1549         findJoinNode();
  1550         bool change = findLeavingArc();
  1551         changeFlow(change);
  1552         if (change) {
  1553           updateTreeStructure();
  1554           updatePotential();
  1555         }
  1556       }
  1557 
  1558       // Copy flow values to _flow_map
  1559       if (_plower) {
  1560         for (int i = 0; i != _arc_num; ++i) {
  1561           Arc e = _arc_ref[i];
  1562           _flow_map->set(e, (*_plower)[e] + _flow[i]);
  1563         }
  1564       } else {
  1565         for (int i = 0; i != _arc_num; ++i) {
  1566           _flow_map->set(_arc_ref[i], _flow[i]);
  1567         }
  1568       }
  1569       // Copy potential values to _potential_map
  1570       for (NodeIt n(_graph); n != INVALID; ++n) {
  1571         _potential_map->set(n, _pi[_node_id[n]]);
  1572       }
  1573 
  1574       return true;
  1575     }
  1576 
  1577   }; //class NetworkSimplex
  1578 
  1579   ///@}
  1580 
  1581 } //namespace lemon
  1582 
  1583 #endif //LEMON_NETWORK_SIMPLEX_H