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