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