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