lemon/network_simplex.h
author Alpar Juttner <alpar@cs.elte.hu>
Tue, 28 Apr 2009 14:01:00 +0100
changeset 628 586b65073025
parent 612 0c8e5c688440
parent 613 b1811c363299
child 640 6c408d864fa1
permissions -rw-r--r--
Merge
     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 *
   385                                     std::sqrt(double(_arc_num))),
   386                                 MIN_BLOCK_SIZE );
   387       }
   388 
   389       // Find next entering arc
   390       bool findEnteringArc() {
   391         Cost c, min = 0;
   392         int cnt = _block_size;
   393         int e, min_arc = _next_arc;
   394         for (e = _next_arc; e < _arc_num; ++e) {
   395           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   396           if (c < min) {
   397             min = c;
   398             min_arc = e;
   399           }
   400           if (--cnt == 0) {
   401             if (min < 0) break;
   402             cnt = _block_size;
   403           }
   404         }
   405         if (min == 0 || cnt > 0) {
   406           for (e = 0; e < _next_arc; ++e) {
   407             c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   408             if (c < min) {
   409               min = c;
   410               min_arc = e;
   411             }
   412             if (--cnt == 0) {
   413               if (min < 0) break;
   414               cnt = _block_size;
   415             }
   416           }
   417         }
   418         if (min >= 0) return false;
   419         _in_arc = min_arc;
   420         _next_arc = e;
   421         return true;
   422       }
   423 
   424     }; //class BlockSearchPivotRule
   425 
   426 
   427     // Implementation of the Candidate List pivot rule
   428     class CandidateListPivotRule
   429     {
   430     private:
   431 
   432       // References to the NetworkSimplex class
   433       const IntVector  &_source;
   434       const IntVector  &_target;
   435       const CostVector &_cost;
   436       const IntVector  &_state;
   437       const CostVector &_pi;
   438       int &_in_arc;
   439       int _arc_num;
   440 
   441       // Pivot rule data
   442       IntVector _candidates;
   443       int _list_length, _minor_limit;
   444       int _curr_length, _minor_count;
   445       int _next_arc;
   446 
   447     public:
   448 
   449       /// Constructor
   450       CandidateListPivotRule(NetworkSimplex &ns) :
   451         _source(ns._source), _target(ns._target),
   452         _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   453         _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
   454       {
   455         // The main parameters of the pivot rule
   456         const double LIST_LENGTH_FACTOR = 1.0;
   457         const int MIN_LIST_LENGTH = 10;
   458         const double MINOR_LIMIT_FACTOR = 0.1;
   459         const int MIN_MINOR_LIMIT = 3;
   460 
   461         _list_length = std::max( int(LIST_LENGTH_FACTOR *
   462                                      std::sqrt(double(_arc_num))),
   463                                  MIN_LIST_LENGTH );
   464         _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
   465                                  MIN_MINOR_LIMIT );
   466         _curr_length = _minor_count = 0;
   467         _candidates.resize(_list_length);
   468       }
   469 
   470       /// Find next entering arc
   471       bool findEnteringArc() {
   472         Cost min, c;
   473         int e, min_arc = _next_arc;
   474         if (_curr_length > 0 && _minor_count < _minor_limit) {
   475           // Minor iteration: select the best eligible arc from the
   476           // current candidate list
   477           ++_minor_count;
   478           min = 0;
   479           for (int i = 0; i < _curr_length; ++i) {
   480             e = _candidates[i];
   481             c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   482             if (c < min) {
   483               min = c;
   484               min_arc = e;
   485             }
   486             if (c >= 0) {
   487               _candidates[i--] = _candidates[--_curr_length];
   488             }
   489           }
   490           if (min < 0) {
   491             _in_arc = min_arc;
   492             return true;
   493           }
   494         }
   495 
   496         // Major iteration: build a new candidate list
   497         min = 0;
   498         _curr_length = 0;
   499         for (e = _next_arc; e < _arc_num; ++e) {
   500           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   501           if (c < 0) {
   502             _candidates[_curr_length++] = e;
   503             if (c < min) {
   504               min = c;
   505               min_arc = e;
   506             }
   507             if (_curr_length == _list_length) break;
   508           }
   509         }
   510         if (_curr_length < _list_length) {
   511           for (e = 0; e < _next_arc; ++e) {
   512             c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   513             if (c < 0) {
   514               _candidates[_curr_length++] = e;
   515               if (c < min) {
   516                 min = c;
   517                 min_arc = e;
   518               }
   519               if (_curr_length == _list_length) break;
   520             }
   521           }
   522         }
   523         if (_curr_length == 0) return false;
   524         _minor_count = 1;
   525         _in_arc = min_arc;
   526         _next_arc = e;
   527         return true;
   528       }
   529 
   530     }; //class CandidateListPivotRule
   531 
   532 
   533     // Implementation of the Altering Candidate List pivot rule
   534     class AlteringListPivotRule
   535     {
   536     private:
   537 
   538       // References to the NetworkSimplex class
   539       const IntVector  &_source;
   540       const IntVector  &_target;
   541       const CostVector &_cost;
   542       const IntVector  &_state;
   543       const CostVector &_pi;
   544       int &_in_arc;
   545       int _arc_num;
   546 
   547       // Pivot rule data
   548       int _block_size, _head_length, _curr_length;
   549       int _next_arc;
   550       IntVector _candidates;
   551       CostVector _cand_cost;
   552 
   553       // Functor class to compare arcs during sort of the candidate list
   554       class SortFunc
   555       {
   556       private:
   557         const CostVector &_map;
   558       public:
   559         SortFunc(const CostVector &map) : _map(map) {}
   560         bool operator()(int left, int right) {
   561           return _map[left] > _map[right];
   562         }
   563       };
   564 
   565       SortFunc _sort_func;
   566 
   567     public:
   568 
   569       // Constructor
   570       AlteringListPivotRule(NetworkSimplex &ns) :
   571         _source(ns._source), _target(ns._target),
   572         _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   573         _in_arc(ns.in_arc), _arc_num(ns._arc_num),
   574         _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
   575       {
   576         // The main parameters of the pivot rule
   577         const double BLOCK_SIZE_FACTOR = 1.5;
   578         const int MIN_BLOCK_SIZE = 10;
   579         const double HEAD_LENGTH_FACTOR = 0.1;
   580         const int MIN_HEAD_LENGTH = 3;
   581 
   582         _block_size = std::max( int(BLOCK_SIZE_FACTOR *
   583                                     std::sqrt(double(_arc_num))),
   584                                 MIN_BLOCK_SIZE );
   585         _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
   586                                  MIN_HEAD_LENGTH );
   587         _candidates.resize(_head_length + _block_size);
   588         _curr_length = 0;
   589       }
   590 
   591       // Find next entering arc
   592       bool findEnteringArc() {
   593         // Check the current candidate list
   594         int e;
   595         for (int i = 0; i < _curr_length; ++i) {
   596           e = _candidates[i];
   597           _cand_cost[e] = _state[e] *
   598             (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   599           if (_cand_cost[e] >= 0) {
   600             _candidates[i--] = _candidates[--_curr_length];
   601           }
   602         }
   603 
   604         // Extend the list
   605         int cnt = _block_size;
   606         int last_arc = 0;
   607         int limit = _head_length;
   608 
   609         for (int e = _next_arc; e < _arc_num; ++e) {
   610           _cand_cost[e] = _state[e] *
   611             (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   612           if (_cand_cost[e] < 0) {
   613             _candidates[_curr_length++] = e;
   614             last_arc = e;
   615           }
   616           if (--cnt == 0) {
   617             if (_curr_length > limit) break;
   618             limit = 0;
   619             cnt = _block_size;
   620           }
   621         }
   622         if (_curr_length <= limit) {
   623           for (int e = 0; e < _next_arc; ++e) {
   624             _cand_cost[e] = _state[e] *
   625               (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   626             if (_cand_cost[e] < 0) {
   627               _candidates[_curr_length++] = e;
   628               last_arc = e;
   629             }
   630             if (--cnt == 0) {
   631               if (_curr_length > limit) break;
   632               limit = 0;
   633               cnt = _block_size;
   634             }
   635           }
   636         }
   637         if (_curr_length == 0) return false;
   638         _next_arc = last_arc + 1;
   639 
   640         // Make heap of the candidate list (approximating a partial sort)
   641         make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
   642                    _sort_func );
   643 
   644         // Pop the first element of the heap
   645         _in_arc = _candidates[0];
   646         pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
   647                   _sort_func );
   648         _curr_length = std::min(_head_length, _curr_length - 1);
   649         return true;
   650       }
   651 
   652     }; //class AlteringListPivotRule
   653 
   654   public:
   655 
   656     /// \brief Constructor.
   657     ///
   658     /// The constructor of the class.
   659     ///
   660     /// \param graph The digraph the algorithm runs on.
   661     NetworkSimplex(const GR& graph) :
   662       _graph(graph),
   663       _plower(NULL), _pupper(NULL), _pcost(NULL),
   664       _psupply(NULL), _pstsup(false), _ptype(GEQ),
   665       _flow_map(NULL), _potential_map(NULL),
   666       _local_flow(false), _local_potential(false),
   667       _node_id(graph)
   668     {
   669       LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
   670                    std::numeric_limits<Flow>::is_signed,
   671         "The flow type of NetworkSimplex must be signed integer");
   672       LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
   673                    std::numeric_limits<Cost>::is_signed,
   674         "The cost type of NetworkSimplex must be signed integer");
   675     }
   676 
   677     /// Destructor.
   678     ~NetworkSimplex() {
   679       if (_local_flow) delete _flow_map;
   680       if (_local_potential) delete _potential_map;
   681     }
   682 
   683     /// \name Parameters
   684     /// The parameters of the algorithm can be specified using these
   685     /// functions.
   686 
   687     /// @{
   688 
   689     /// \brief Set the lower bounds on the arcs.
   690     ///
   691     /// This function sets the lower bounds on the arcs.
   692     /// If neither this function nor \ref boundMaps() is used before
   693     /// calling \ref run(), the lower bounds will be set to zero
   694     /// on all arcs.
   695     ///
   696     /// \param map An arc map storing the lower bounds.
   697     /// Its \c Value type must be convertible to the \c Flow type
   698     /// of the algorithm.
   699     ///
   700     /// \return <tt>(*this)</tt>
   701     template <typename LOWER>
   702     NetworkSimplex& lowerMap(const LOWER& map) {
   703       delete _plower;
   704       _plower = new FlowArcMap(_graph);
   705       for (ArcIt a(_graph); a != INVALID; ++a) {
   706         (*_plower)[a] = map[a];
   707       }
   708       return *this;
   709     }
   710 
   711     /// \brief Set the upper bounds (capacities) on the arcs.
   712     ///
   713     /// This function sets the upper bounds (capacities) on the arcs.
   714     /// If none of the functions \ref upperMap(), \ref capacityMap()
   715     /// and \ref boundMaps() is used before calling \ref run(),
   716     /// the upper bounds (capacities) will be set to
   717     /// \c std::numeric_limits<Flow>::max() on all arcs.
   718     ///
   719     /// \param map An arc map storing the upper bounds.
   720     /// Its \c Value type must be convertible to the \c Flow type
   721     /// of the algorithm.
   722     ///
   723     /// \return <tt>(*this)</tt>
   724     template<typename UPPER>
   725     NetworkSimplex& upperMap(const UPPER& map) {
   726       delete _pupper;
   727       _pupper = new FlowArcMap(_graph);
   728       for (ArcIt a(_graph); a != INVALID; ++a) {
   729         (*_pupper)[a] = map[a];
   730       }
   731       return *this;
   732     }
   733 
   734     /// \brief Set the upper bounds (capacities) on the arcs.
   735     ///
   736     /// This function sets the upper bounds (capacities) on the arcs.
   737     /// It is just an alias for \ref upperMap().
   738     ///
   739     /// \return <tt>(*this)</tt>
   740     template<typename CAP>
   741     NetworkSimplex& capacityMap(const CAP& map) {
   742       return upperMap(map);
   743     }
   744 
   745     /// \brief Set the lower and upper bounds on the arcs.
   746     ///
   747     /// This function sets the lower and upper bounds on the arcs.
   748     /// If neither this function nor \ref lowerMap() is used before
   749     /// calling \ref run(), the lower bounds will be set to zero
   750     /// on all arcs.
   751     /// If none of the functions \ref upperMap(), \ref capacityMap()
   752     /// and \ref boundMaps() is used before calling \ref run(),
   753     /// the upper bounds (capacities) will be set to
   754     /// \c std::numeric_limits<Flow>::max() on all arcs.
   755     ///
   756     /// \param lower An arc map storing the lower bounds.
   757     /// \param upper An arc map storing the upper bounds.
   758     ///
   759     /// The \c Value type of the maps must be convertible to the
   760     /// \c Flow type of the algorithm.
   761     ///
   762     /// \note This function is just a shortcut of calling \ref lowerMap()
   763     /// and \ref upperMap() separately.
   764     ///
   765     /// \return <tt>(*this)</tt>
   766     template <typename LOWER, typename UPPER>
   767     NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
   768       return lowerMap(lower).upperMap(upper);
   769     }
   770 
   771     /// \brief Set the costs of the arcs.
   772     ///
   773     /// This function sets the costs of the arcs.
   774     /// If it is not used before calling \ref run(), the costs
   775     /// will be set to \c 1 on all arcs.
   776     ///
   777     /// \param map An arc map storing the costs.
   778     /// Its \c Value type must be convertible to the \c Cost type
   779     /// of the algorithm.
   780     ///
   781     /// \return <tt>(*this)</tt>
   782     template<typename COST>
   783     NetworkSimplex& costMap(const COST& map) {
   784       delete _pcost;
   785       _pcost = new CostArcMap(_graph);
   786       for (ArcIt a(_graph); a != INVALID; ++a) {
   787         (*_pcost)[a] = map[a];
   788       }
   789       return *this;
   790     }
   791 
   792     /// \brief Set the supply values of the nodes.
   793     ///
   794     /// This function sets the supply values of the nodes.
   795     /// If neither this function nor \ref stSupply() is used before
   796     /// calling \ref run(), the supply of each node will be set to zero.
   797     /// (It makes sense only if non-zero lower bounds are given.)
   798     ///
   799     /// \param map A node map storing the supply values.
   800     /// Its \c Value type must be convertible to the \c Flow type
   801     /// of the algorithm.
   802     ///
   803     /// \return <tt>(*this)</tt>
   804     template<typename SUP>
   805     NetworkSimplex& supplyMap(const SUP& map) {
   806       delete _psupply;
   807       _pstsup = false;
   808       _psupply = new FlowNodeMap(_graph);
   809       for (NodeIt n(_graph); n != INVALID; ++n) {
   810         (*_psupply)[n] = map[n];
   811       }
   812       return *this;
   813     }
   814 
   815     /// \brief Set single source and target nodes and a supply value.
   816     ///
   817     /// This function sets a single source node and a single target node
   818     /// and the required flow value.
   819     /// If neither this function nor \ref supplyMap() is used before
   820     /// calling \ref run(), the supply of each node will be set to zero.
   821     /// (It makes sense only if non-zero lower bounds are given.)
   822     ///
   823     /// \param s The source node.
   824     /// \param t The target node.
   825     /// \param k The required amount of flow from node \c s to node \c t
   826     /// (i.e. the supply of \c s and the demand of \c t).
   827     ///
   828     /// \return <tt>(*this)</tt>
   829     NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) {
   830       delete _psupply;
   831       _psupply = NULL;
   832       _pstsup = true;
   833       _psource = s;
   834       _ptarget = t;
   835       _pstflow = k;
   836       return *this;
   837     }
   838     
   839     /// \brief Set the problem type.
   840     ///
   841     /// This function sets the problem type for the algorithm.
   842     /// If it is not used before calling \ref run(), the \ref GEQ problem
   843     /// type will be used.
   844     ///
   845     /// For more information see \ref ProblemType.
   846     ///
   847     /// \return <tt>(*this)</tt>
   848     NetworkSimplex& problemType(ProblemType problem_type) {
   849       _ptype = problem_type;
   850       return *this;
   851     }
   852 
   853     /// \brief Set the flow map.
   854     ///
   855     /// This function sets the flow map.
   856     /// If it is not used before calling \ref run(), an instance will
   857     /// be allocated automatically. The destructor deallocates this
   858     /// automatically allocated map, of course.
   859     ///
   860     /// \return <tt>(*this)</tt>
   861     NetworkSimplex& flowMap(FlowMap& map) {
   862       if (_local_flow) {
   863         delete _flow_map;
   864         _local_flow = false;
   865       }
   866       _flow_map = &map;
   867       return *this;
   868     }
   869 
   870     /// \brief Set the potential map.
   871     ///
   872     /// This function sets the potential map, which is used for storing
   873     /// the dual solution.
   874     /// If it is not used before calling \ref run(), an instance will
   875     /// be allocated automatically. The destructor deallocates this
   876     /// automatically allocated map, of course.
   877     ///
   878     /// \return <tt>(*this)</tt>
   879     NetworkSimplex& potentialMap(PotentialMap& map) {
   880       if (_local_potential) {
   881         delete _potential_map;
   882         _local_potential = false;
   883       }
   884       _potential_map = &map;
   885       return *this;
   886     }
   887     
   888     /// @}
   889 
   890     /// \name Execution Control
   891     /// The algorithm can be executed using \ref run().
   892 
   893     /// @{
   894 
   895     /// \brief Run the algorithm.
   896     ///
   897     /// This function runs the algorithm.
   898     /// The paramters can be specified using functions \ref lowerMap(),
   899     /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
   900     /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 
   901     /// \ref problemType(), \ref flowMap() and \ref potentialMap().
   902     /// For example,
   903     /// \code
   904     ///   NetworkSimplex<ListDigraph> ns(graph);
   905     ///   ns.boundMaps(lower, upper).costMap(cost)
   906     ///     .supplyMap(sup).run();
   907     /// \endcode
   908     ///
   909     /// This function can be called more than once. All the parameters
   910     /// that have been given are kept for the next call, unless
   911     /// \ref reset() is called, thus only the modified parameters
   912     /// have to be set again. See \ref reset() for examples.
   913     ///
   914     /// \param pivot_rule The pivot rule that will be used during the
   915     /// algorithm. For more information see \ref PivotRule.
   916     ///
   917     /// \return \c true if a feasible flow can be found.
   918     bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
   919       return init() && start(pivot_rule);
   920     }
   921 
   922     /// \brief Reset all the parameters that have been given before.
   923     ///
   924     /// This function resets all the paramaters that have been given
   925     /// before using functions \ref lowerMap(), \ref upperMap(),
   926     /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
   927     /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 
   928     /// \ref flowMap() and \ref potentialMap().
   929     ///
   930     /// It is useful for multiple run() calls. If this function is not
   931     /// used, all the parameters given before are kept for the next
   932     /// \ref run() call.
   933     ///
   934     /// For example,
   935     /// \code
   936     ///   NetworkSimplex<ListDigraph> ns(graph);
   937     ///
   938     ///   // First run
   939     ///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)
   940     ///     .supplyMap(sup).run();
   941     ///
   942     ///   // Run again with modified cost map (reset() is not called,
   943     ///   // so only the cost map have to be set again)
   944     ///   cost[e] += 100;
   945     ///   ns.costMap(cost).run();
   946     ///
   947     ///   // Run again from scratch using reset()
   948     ///   // (the lower bounds will be set to zero on all arcs)
   949     ///   ns.reset();
   950     ///   ns.capacityMap(cap).costMap(cost)
   951     ///     .supplyMap(sup).run();
   952     /// \endcode
   953     ///
   954     /// \return <tt>(*this)</tt>
   955     NetworkSimplex& reset() {
   956       delete _plower;
   957       delete _pupper;
   958       delete _pcost;
   959       delete _psupply;
   960       _plower = NULL;
   961       _pupper = NULL;
   962       _pcost = NULL;
   963       _psupply = NULL;
   964       _pstsup = false;
   965       _ptype = GEQ;
   966       if (_local_flow) delete _flow_map;
   967       if (_local_potential) delete _potential_map;
   968       _flow_map = NULL;
   969       _potential_map = NULL;
   970       _local_flow = false;
   971       _local_potential = false;
   972 
   973       return *this;
   974     }
   975 
   976     /// @}
   977 
   978     /// \name Query Functions
   979     /// The results of the algorithm can be obtained using these
   980     /// functions.\n
   981     /// The \ref run() function must be called before using them.
   982 
   983     /// @{
   984 
   985     /// \brief Return the total cost of the found flow.
   986     ///
   987     /// This function returns the total cost of the found flow.
   988     /// The complexity of the function is O(e).
   989     ///
   990     /// \note The return type of the function can be specified as a
   991     /// template parameter. For example,
   992     /// \code
   993     ///   ns.totalCost<double>();
   994     /// \endcode
   995     /// It is useful if the total cost cannot be stored in the \c Cost
   996     /// type of the algorithm, which is the default return type of the
   997     /// function.
   998     ///
   999     /// \pre \ref run() must be called before using this function.
  1000     template <typename Num>
  1001     Num totalCost() const {
  1002       Num c = 0;
  1003       if (_pcost) {
  1004         for (ArcIt e(_graph); e != INVALID; ++e)
  1005           c += (*_flow_map)[e] * (*_pcost)[e];
  1006       } else {
  1007         for (ArcIt e(_graph); e != INVALID; ++e)
  1008           c += (*_flow_map)[e];
  1009       }
  1010       return c;
  1011     }
  1012 
  1013 #ifndef DOXYGEN
  1014     Cost totalCost() const {
  1015       return totalCost<Cost>();
  1016     }
  1017 #endif
  1018 
  1019     /// \brief Return the flow on the given arc.
  1020     ///
  1021     /// This function returns the flow on the given arc.
  1022     ///
  1023     /// \pre \ref run() must be called before using this function.
  1024     Flow flow(const Arc& a) const {
  1025       return (*_flow_map)[a];
  1026     }
  1027 
  1028     /// \brief Return a const reference to the flow map.
  1029     ///
  1030     /// This function returns a const reference to an arc map storing
  1031     /// the found flow.
  1032     ///
  1033     /// \pre \ref run() must be called before using this function.
  1034     const FlowMap& flowMap() const {
  1035       return *_flow_map;
  1036     }
  1037 
  1038     /// \brief Return the potential (dual value) of the given node.
  1039     ///
  1040     /// This function returns the potential (dual value) of the
  1041     /// given node.
  1042     ///
  1043     /// \pre \ref run() must be called before using this function.
  1044     Cost potential(const Node& n) const {
  1045       return (*_potential_map)[n];
  1046     }
  1047 
  1048     /// \brief Return a const reference to the potential map
  1049     /// (the dual solution).
  1050     ///
  1051     /// This function returns a const reference to a node map storing
  1052     /// the found potentials, which form the dual solution of the
  1053     /// \ref min_cost_flow "minimum cost flow" problem.
  1054     ///
  1055     /// \pre \ref run() must be called before using this function.
  1056     const PotentialMap& potentialMap() const {
  1057       return *_potential_map;
  1058     }
  1059 
  1060     /// @}
  1061 
  1062   private:
  1063 
  1064     // Initialize internal data structures
  1065     bool init() {
  1066       // Initialize result maps
  1067       if (!_flow_map) {
  1068         _flow_map = new FlowMap(_graph);
  1069         _local_flow = true;
  1070       }
  1071       if (!_potential_map) {
  1072         _potential_map = new PotentialMap(_graph);
  1073         _local_potential = true;
  1074       }
  1075 
  1076       // Initialize vectors
  1077       _node_num = countNodes(_graph);
  1078       _arc_num = countArcs(_graph);
  1079       int all_node_num = _node_num + 1;
  1080       int all_arc_num = _arc_num + _node_num;
  1081       if (_node_num == 0) return false;
  1082 
  1083       _arc_ref.resize(_arc_num);
  1084       _source.resize(all_arc_num);
  1085       _target.resize(all_arc_num);
  1086 
  1087       _cap.resize(all_arc_num);
  1088       _cost.resize(all_arc_num);
  1089       _supply.resize(all_node_num);
  1090       _flow.resize(all_arc_num);
  1091       _pi.resize(all_node_num);
  1092 
  1093       _parent.resize(all_node_num);
  1094       _pred.resize(all_node_num);
  1095       _forward.resize(all_node_num);
  1096       _thread.resize(all_node_num);
  1097       _rev_thread.resize(all_node_num);
  1098       _succ_num.resize(all_node_num);
  1099       _last_succ.resize(all_node_num);
  1100       _state.resize(all_arc_num);
  1101 
  1102       // Initialize node related data
  1103       bool valid_supply = true;
  1104       Flow sum_supply = 0;
  1105       if (!_pstsup && !_psupply) {
  1106         _pstsup = true;
  1107         _psource = _ptarget = NodeIt(_graph);
  1108         _pstflow = 0;
  1109       }
  1110       if (_psupply) {
  1111         int i = 0;
  1112         for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
  1113           _node_id[n] = i;
  1114           _supply[i] = (*_psupply)[n];
  1115           sum_supply += _supply[i];
  1116         }
  1117         valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
  1118                        (_ptype == LEQ && sum_supply >= 0);
  1119       } else {
  1120         int i = 0;
  1121         for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
  1122           _node_id[n] = i;
  1123           _supply[i] = 0;
  1124         }
  1125         _supply[_node_id[_psource]] =  _pstflow;
  1126         _supply[_node_id[_ptarget]] = -_pstflow;
  1127       }
  1128       if (!valid_supply) return false;
  1129 
  1130       // Infinite capacity value
  1131       Flow inf_cap =
  1132         std::numeric_limits<Flow>::has_infinity ?
  1133         std::numeric_limits<Flow>::infinity() :
  1134         std::numeric_limits<Flow>::max();
  1135 
  1136       // Initialize artifical cost
  1137       Cost art_cost;
  1138       if (std::numeric_limits<Cost>::is_exact) {
  1139         art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
  1140       } else {
  1141         art_cost = std::numeric_limits<Cost>::min();
  1142         for (int i = 0; i != _arc_num; ++i) {
  1143           if (_cost[i] > art_cost) art_cost = _cost[i];
  1144         }
  1145         art_cost = (art_cost + 1) * _node_num;
  1146       }
  1147 
  1148       // Run Circulation to check if a feasible solution exists
  1149       typedef ConstMap<Arc, Flow> ConstArcMap;
  1150       ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap);
  1151       FlowNodeMap *csup = NULL;
  1152       bool local_csup = false;
  1153       if (_psupply) {
  1154         csup = _psupply;
  1155       } else {
  1156         csup = new FlowNodeMap(_graph, 0);
  1157         (*csup)[_psource] =  _pstflow;
  1158         (*csup)[_ptarget] = -_pstflow;
  1159         local_csup = true;
  1160       }
  1161       bool circ_result = false;
  1162       if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
  1163         // GEQ problem type
  1164         if (_plower) {
  1165           if (_pupper) {
  1166             Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
  1167               circ(_graph, *_plower, *_pupper, *csup);
  1168             circ_result = circ.run();
  1169           } else {
  1170             Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
  1171               circ(_graph, *_plower, inf_arc_map, *csup);
  1172             circ_result = circ.run();
  1173           }
  1174         } else {
  1175           if (_pupper) {
  1176             Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
  1177               circ(_graph, zero_arc_map, *_pupper, *csup);
  1178             circ_result = circ.run();
  1179           } else {
  1180             Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
  1181               circ(_graph, zero_arc_map, inf_arc_map, *csup);
  1182             circ_result = circ.run();
  1183           }
  1184         }
  1185       } else {
  1186         // LEQ problem type
  1187         typedef ReverseDigraph<const GR> RevGraph;
  1188         typedef NegMap<FlowNodeMap> NegNodeMap;
  1189         RevGraph rgraph(_graph);
  1190         NegNodeMap neg_csup(*csup);
  1191         if (_plower) {
  1192           if (_pupper) {
  1193             Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
  1194               circ(rgraph, *_plower, *_pupper, neg_csup);
  1195             circ_result = circ.run();
  1196           } else {
  1197             Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
  1198               circ(rgraph, *_plower, inf_arc_map, neg_csup);
  1199             circ_result = circ.run();
  1200           }
  1201         } else {
  1202           if (_pupper) {
  1203             Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
  1204               circ(rgraph, zero_arc_map, *_pupper, neg_csup);
  1205             circ_result = circ.run();
  1206           } else {
  1207             Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
  1208               circ(rgraph, zero_arc_map, inf_arc_map, neg_csup);
  1209             circ_result = circ.run();
  1210           }
  1211         }
  1212       }
  1213       if (local_csup) delete csup;
  1214       if (!circ_result) return false;
  1215 
  1216       // Set data for the artificial root node
  1217       _root = _node_num;
  1218       _parent[_root] = -1;
  1219       _pred[_root] = -1;
  1220       _thread[_root] = 0;
  1221       _rev_thread[0] = _root;
  1222       _succ_num[_root] = all_node_num;
  1223       _last_succ[_root] = _root - 1;
  1224       _supply[_root] = -sum_supply;
  1225       if (sum_supply < 0) {
  1226         _pi[_root] = -art_cost;
  1227       } else {
  1228         _pi[_root] = art_cost;
  1229       }
  1230 
  1231       // Store the arcs in a mixed order
  1232       int k = std::max(int(std::sqrt(double(_arc_num))), 10);
  1233       int i = 0;
  1234       for (ArcIt e(_graph); e != INVALID; ++e) {
  1235         _arc_ref[i] = e;
  1236         if ((i += k) >= _arc_num) i = (i % k) + 1;
  1237       }
  1238 
  1239       // Initialize arc maps
  1240       if (_pupper && _pcost) {
  1241         for (int i = 0; i != _arc_num; ++i) {
  1242           Arc e = _arc_ref[i];
  1243           _source[i] = _node_id[_graph.source(e)];
  1244           _target[i] = _node_id[_graph.target(e)];
  1245           _cap[i] = (*_pupper)[e];
  1246           _cost[i] = (*_pcost)[e];
  1247           _flow[i] = 0;
  1248           _state[i] = STATE_LOWER;
  1249         }
  1250       } else {
  1251         for (int i = 0; i != _arc_num; ++i) {
  1252           Arc e = _arc_ref[i];
  1253           _source[i] = _node_id[_graph.source(e)];
  1254           _target[i] = _node_id[_graph.target(e)];
  1255           _flow[i] = 0;
  1256           _state[i] = STATE_LOWER;
  1257         }
  1258         if (_pupper) {
  1259           for (int i = 0; i != _arc_num; ++i)
  1260             _cap[i] = (*_pupper)[_arc_ref[i]];
  1261         } else {
  1262           for (int i = 0; i != _arc_num; ++i)
  1263             _cap[i] = inf_cap;
  1264         }
  1265         if (_pcost) {
  1266           for (int i = 0; i != _arc_num; ++i)
  1267             _cost[i] = (*_pcost)[_arc_ref[i]];
  1268         } else {
  1269           for (int i = 0; i != _arc_num; ++i)
  1270             _cost[i] = 1;
  1271         }
  1272       }
  1273       
  1274       // Remove non-zero lower bounds
  1275       if (_plower) {
  1276         for (int i = 0; i != _arc_num; ++i) {
  1277           Flow c = (*_plower)[_arc_ref[i]];
  1278           if (c != 0) {
  1279             _cap[i] -= c;
  1280             _supply[_source[i]] -= c;
  1281             _supply[_target[i]] += c;
  1282           }
  1283         }
  1284       }
  1285 
  1286       // Add artificial arcs and initialize the spanning tree data structure
  1287       for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
  1288         _thread[u] = u + 1;
  1289         _rev_thread[u + 1] = u;
  1290         _succ_num[u] = 1;
  1291         _last_succ[u] = u;
  1292         _parent[u] = _root;
  1293         _pred[u] = e;
  1294         _cost[e] = art_cost;
  1295         _cap[e] = inf_cap;
  1296         _state[e] = STATE_TREE;
  1297         if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
  1298           _flow[e] = _supply[u];
  1299           _forward[u] = true;
  1300           _pi[u] = -art_cost + _pi[_root];
  1301         } else {
  1302           _flow[e] = -_supply[u];
  1303           _forward[u] = false;
  1304           _pi[u] = art_cost + _pi[_root];
  1305         }
  1306       }
  1307 
  1308       return true;
  1309     }
  1310 
  1311     // Find the join node
  1312     void findJoinNode() {
  1313       int u = _source[in_arc];
  1314       int v = _target[in_arc];
  1315       while (u != v) {
  1316         if (_succ_num[u] < _succ_num[v]) {
  1317           u = _parent[u];
  1318         } else {
  1319           v = _parent[v];
  1320         }
  1321       }
  1322       join = u;
  1323     }
  1324 
  1325     // Find the leaving arc of the cycle and returns true if the
  1326     // leaving arc is not the same as the entering arc
  1327     bool findLeavingArc() {
  1328       // Initialize first and second nodes according to the direction
  1329       // of the cycle
  1330       if (_state[in_arc] == STATE_LOWER) {
  1331         first  = _source[in_arc];
  1332         second = _target[in_arc];
  1333       } else {
  1334         first  = _target[in_arc];
  1335         second = _source[in_arc];
  1336       }
  1337       delta = _cap[in_arc];
  1338       int result = 0;
  1339       Flow d;
  1340       int e;
  1341 
  1342       // Search the cycle along the path form the first node to the root
  1343       for (int u = first; u != join; u = _parent[u]) {
  1344         e = _pred[u];
  1345         d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
  1346         if (d < delta) {
  1347           delta = d;
  1348           u_out = u;
  1349           result = 1;
  1350         }
  1351       }
  1352       // Search the cycle along the path form the second node to the root
  1353       for (int u = second; u != join; u = _parent[u]) {
  1354         e = _pred[u];
  1355         d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
  1356         if (d <= delta) {
  1357           delta = d;
  1358           u_out = u;
  1359           result = 2;
  1360         }
  1361       }
  1362 
  1363       if (result == 1) {
  1364         u_in = first;
  1365         v_in = second;
  1366       } else {
  1367         u_in = second;
  1368         v_in = first;
  1369       }
  1370       return result != 0;
  1371     }
  1372 
  1373     // Change _flow and _state vectors
  1374     void changeFlow(bool change) {
  1375       // Augment along the cycle
  1376       if (delta > 0) {
  1377         Flow val = _state[in_arc] * delta;
  1378         _flow[in_arc] += val;
  1379         for (int u = _source[in_arc]; u != join; u = _parent[u]) {
  1380           _flow[_pred[u]] += _forward[u] ? -val : val;
  1381         }
  1382         for (int u = _target[in_arc]; u != join; u = _parent[u]) {
  1383           _flow[_pred[u]] += _forward[u] ? val : -val;
  1384         }
  1385       }
  1386       // Update the state of the entering and leaving arcs
  1387       if (change) {
  1388         _state[in_arc] = STATE_TREE;
  1389         _state[_pred[u_out]] =
  1390           (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
  1391       } else {
  1392         _state[in_arc] = -_state[in_arc];
  1393       }
  1394     }
  1395 
  1396     // Update the tree structure
  1397     void updateTreeStructure() {
  1398       int u, w;
  1399       int old_rev_thread = _rev_thread[u_out];
  1400       int old_succ_num = _succ_num[u_out];
  1401       int old_last_succ = _last_succ[u_out];
  1402       v_out = _parent[u_out];
  1403 
  1404       u = _last_succ[u_in];  // the last successor of u_in
  1405       right = _thread[u];    // the node after it
  1406 
  1407       // Handle the case when old_rev_thread equals to v_in
  1408       // (it also means that join and v_out coincide)
  1409       if (old_rev_thread == v_in) {
  1410         last = _thread[_last_succ[u_out]];
  1411       } else {
  1412         last = _thread[v_in];
  1413       }
  1414 
  1415       // Update _thread and _parent along the stem nodes (i.e. the nodes
  1416       // between u_in and u_out, whose parent have to be changed)
  1417       _thread[v_in] = stem = u_in;
  1418       _dirty_revs.clear();
  1419       _dirty_revs.push_back(v_in);
  1420       par_stem = v_in;
  1421       while (stem != u_out) {
  1422         // Insert the next stem node into the thread list
  1423         new_stem = _parent[stem];
  1424         _thread[u] = new_stem;
  1425         _dirty_revs.push_back(u);
  1426 
  1427         // Remove the subtree of stem from the thread list
  1428         w = _rev_thread[stem];
  1429         _thread[w] = right;
  1430         _rev_thread[right] = w;
  1431 
  1432         // Change the parent node and shift stem nodes
  1433         _parent[stem] = par_stem;
  1434         par_stem = stem;
  1435         stem = new_stem;
  1436 
  1437         // Update u and right
  1438         u = _last_succ[stem] == _last_succ[par_stem] ?
  1439           _rev_thread[par_stem] : _last_succ[stem];
  1440         right = _thread[u];
  1441       }
  1442       _parent[u_out] = par_stem;
  1443       _thread[u] = last;
  1444       _rev_thread[last] = u;
  1445       _last_succ[u_out] = u;
  1446 
  1447       // Remove the subtree of u_out from the thread list except for
  1448       // the case when old_rev_thread equals to v_in
  1449       // (it also means that join and v_out coincide)
  1450       if (old_rev_thread != v_in) {
  1451         _thread[old_rev_thread] = right;
  1452         _rev_thread[right] = old_rev_thread;
  1453       }
  1454 
  1455       // Update _rev_thread using the new _thread values
  1456       for (int i = 0; i < int(_dirty_revs.size()); ++i) {
  1457         u = _dirty_revs[i];
  1458         _rev_thread[_thread[u]] = u;
  1459       }
  1460 
  1461       // Update _pred, _forward, _last_succ and _succ_num for the
  1462       // stem nodes from u_out to u_in
  1463       int tmp_sc = 0, tmp_ls = _last_succ[u_out];
  1464       u = u_out;
  1465       while (u != u_in) {
  1466         w = _parent[u];
  1467         _pred[u] = _pred[w];
  1468         _forward[u] = !_forward[w];
  1469         tmp_sc += _succ_num[u] - _succ_num[w];
  1470         _succ_num[u] = tmp_sc;
  1471         _last_succ[w] = tmp_ls;
  1472         u = w;
  1473       }
  1474       _pred[u_in] = in_arc;
  1475       _forward[u_in] = (u_in == _source[in_arc]);
  1476       _succ_num[u_in] = old_succ_num;
  1477 
  1478       // Set limits for updating _last_succ form v_in and v_out
  1479       // towards the root
  1480       int up_limit_in = -1;
  1481       int up_limit_out = -1;
  1482       if (_last_succ[join] == v_in) {
  1483         up_limit_out = join;
  1484       } else {
  1485         up_limit_in = join;
  1486       }
  1487 
  1488       // Update _last_succ from v_in towards the root
  1489       for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
  1490            u = _parent[u]) {
  1491         _last_succ[u] = _last_succ[u_out];
  1492       }
  1493       // Update _last_succ from v_out towards the root
  1494       if (join != old_rev_thread && v_in != old_rev_thread) {
  1495         for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
  1496              u = _parent[u]) {
  1497           _last_succ[u] = old_rev_thread;
  1498         }
  1499       } else {
  1500         for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
  1501              u = _parent[u]) {
  1502           _last_succ[u] = _last_succ[u_out];
  1503         }
  1504       }
  1505 
  1506       // Update _succ_num from v_in to join
  1507       for (u = v_in; u != join; u = _parent[u]) {
  1508         _succ_num[u] += old_succ_num;
  1509       }
  1510       // Update _succ_num from v_out to join
  1511       for (u = v_out; u != join; u = _parent[u]) {
  1512         _succ_num[u] -= old_succ_num;
  1513       }
  1514     }
  1515 
  1516     // Update potentials
  1517     void updatePotential() {
  1518       Cost sigma = _forward[u_in] ?
  1519         _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
  1520         _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
  1521       // Update potentials in the subtree, which has been moved
  1522       int end = _thread[_last_succ[u_in]];
  1523       for (int u = u_in; u != end; u = _thread[u]) {
  1524         _pi[u] += sigma;
  1525       }
  1526     }
  1527 
  1528     // Execute the algorithm
  1529     bool start(PivotRule pivot_rule) {
  1530       // Select the pivot rule implementation
  1531       switch (pivot_rule) {
  1532         case FIRST_ELIGIBLE:
  1533           return start<FirstEligiblePivotRule>();
  1534         case BEST_ELIGIBLE:
  1535           return start<BestEligiblePivotRule>();
  1536         case BLOCK_SEARCH:
  1537           return start<BlockSearchPivotRule>();
  1538         case CANDIDATE_LIST:
  1539           return start<CandidateListPivotRule>();
  1540         case ALTERING_LIST:
  1541           return start<AlteringListPivotRule>();
  1542       }
  1543       return false;
  1544     }
  1545 
  1546     template <typename PivotRuleImpl>
  1547     bool start() {
  1548       PivotRuleImpl pivot(*this);
  1549 
  1550       // Execute the Network Simplex algorithm
  1551       while (pivot.findEnteringArc()) {
  1552         findJoinNode();
  1553         bool change = findLeavingArc();
  1554         changeFlow(change);
  1555         if (change) {
  1556           updateTreeStructure();
  1557           updatePotential();
  1558         }
  1559       }
  1560 
  1561       // Copy flow values to _flow_map
  1562       if (_plower) {
  1563         for (int i = 0; i != _arc_num; ++i) {
  1564           Arc e = _arc_ref[i];
  1565           _flow_map->set(e, (*_plower)[e] + _flow[i]);
  1566         }
  1567       } else {
  1568         for (int i = 0; i != _arc_num; ++i) {
  1569           _flow_map->set(_arc_ref[i], _flow[i]);
  1570         }
  1571       }
  1572       // Copy potential values to _potential_map
  1573       for (NodeIt n(_graph); n != INVALID; ++n) {
  1574         _potential_map->set(n, _pi[_node_id[n]]);
  1575       }
  1576 
  1577       return true;
  1578     }
  1579 
  1580   }; //class NetworkSimplex
  1581 
  1582   ///@}
  1583 
  1584 } //namespace lemon
  1585 
  1586 #endif //LEMON_NETWORK_SIMPLEX_H