lemon/network_simplex.h
 author Alpar Juttner Wed, 29 Jul 2020 14:56:10 +0200 changeset 1433 a278d16bd2d0 parent 1298 a78e5b779b69 parent 1317 b40c2bbb8da5 permissions -rw-r--r--
Fix clang compilation issue (#634)
     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-2013

     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   /// \cite amo93networkflows, \cite dantzig63linearprog,

    45   /// \cite 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 _has_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       _has_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       _has_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 && _node_num > 1) {

   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(m).

   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       // Check lower and upper bounds

  1071       LEMON_DEBUG(checkBoundMaps(),

  1072           "Upper bounds must be greater or equal to the lower bounds");

  1073

  1074       // Remove non-zero lower bounds

  1075       if (_has_lower) {

  1076         for (int i = 0; i != _arc_num; ++i) {

  1077           Value c = _lower[i];

  1078           if (c >= 0) {

  1079             _cap[i] = _upper[i] < MAX ? _upper[i] - c : INF;

  1080           } else {

  1081             _cap[i] = _upper[i] < MAX + c ? _upper[i] - c : INF;

  1082           }

  1083           _supply[_source[i]] -= c;

  1084           _supply[_target[i]] += c;

  1085         }

  1086       } else {

  1087         for (int i = 0; i != _arc_num; ++i) {

  1088           _cap[i] = _upper[i];

  1089         }

  1090       }

  1091

  1092       // Initialize artifical cost

  1093       Cost ART_COST;

  1094       if (std::numeric_limits<Cost>::is_exact) {

  1095         ART_COST = std::numeric_limits<Cost>::max() / 2 + 1;

  1096       } else {

  1097         ART_COST = 0;

  1098         for (int i = 0; i != _arc_num; ++i) {

  1099           if (_cost[i] > ART_COST) ART_COST = _cost[i];

  1100         }

  1101         ART_COST = (ART_COST + 1) * _node_num;

  1102       }

  1103

  1104       // Initialize arc maps

  1105       for (int i = 0; i != _arc_num; ++i) {

  1106         _flow[i] = 0;

  1107         _state[i] = STATE_LOWER;

  1108       }

  1109

  1110       // Set data for the artificial root node

  1111       _root = _node_num;

  1112       _parent[_root] = -1;

  1113       _pred[_root] = -1;

  1114       _thread[_root] = 0;

  1115       _rev_thread[0] = _root;

  1116       _succ_num[_root] = _node_num + 1;

  1117       _last_succ[_root] = _root - 1;

  1118       _supply[_root] = -_sum_supply;

  1119       _pi[_root] = 0;

  1120

  1121       // Add artificial arcs and initialize the spanning tree data structure

  1122       if (_sum_supply == 0) {

  1123         // EQ supply constraints

  1124         _search_arc_num = _arc_num;

  1125         _all_arc_num = _arc_num + _node_num;

  1126         for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {

  1127           _parent[u] = _root;

  1128           _pred[u] = e;

  1129           _thread[u] = u + 1;

  1130           _rev_thread[u + 1] = u;

  1131           _succ_num[u] = 1;

  1132           _last_succ[u] = u;

  1133           _cap[e] = INF;

  1134           _state[e] = STATE_TREE;

  1135           if (_supply[u] >= 0) {

  1136             _pred_dir[u] = DIR_UP;

  1137             _pi[u] = 0;

  1138             _source[e] = u;

  1139             _target[e] = _root;

  1140             _flow[e] = _supply[u];

  1141             _cost[e] = 0;

  1142           } else {

  1143             _pred_dir[u] = DIR_DOWN;

  1144             _pi[u] = ART_COST;

  1145             _source[e] = _root;

  1146             _target[e] = u;

  1147             _flow[e] = -_supply[u];

  1148             _cost[e] = ART_COST;

  1149           }

  1150         }

  1151       }

  1152       else if (_sum_supply > 0) {

  1153         // LEQ supply constraints

  1154         _search_arc_num = _arc_num + _node_num;

  1155         int f = _arc_num + _node_num;

  1156         for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {

  1157           _parent[u] = _root;

  1158           _thread[u] = u + 1;

  1159           _rev_thread[u + 1] = u;

  1160           _succ_num[u] = 1;

  1161           _last_succ[u] = u;

  1162           if (_supply[u] >= 0) {

  1163             _pred_dir[u] = DIR_UP;

  1164             _pi[u] = 0;

  1165             _pred[u] = e;

  1166             _source[e] = u;

  1167             _target[e] = _root;

  1168             _cap[e] = INF;

  1169             _flow[e] = _supply[u];

  1170             _cost[e] = 0;

  1171             _state[e] = STATE_TREE;

  1172           } else {

  1173             _pred_dir[u] = DIR_DOWN;

  1174             _pi[u] = ART_COST;

  1175             _pred[u] = f;

  1176             _source[f] = _root;

  1177             _target[f] = u;

  1178             _cap[f] = INF;

  1179             _flow[f] = -_supply[u];

  1180             _cost[f] = ART_COST;

  1181             _state[f] = STATE_TREE;

  1182             _source[e] = u;

  1183             _target[e] = _root;

  1184             _cap[e] = INF;

  1185             _flow[e] = 0;

  1186             _cost[e] = 0;

  1187             _state[e] = STATE_LOWER;

  1188             ++f;

  1189           }

  1190         }

  1191         _all_arc_num = f;

  1192       }

  1193       else {

  1194         // GEQ supply constraints

  1195         _search_arc_num = _arc_num + _node_num;

  1196         int f = _arc_num + _node_num;

  1197         for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {

  1198           _parent[u] = _root;

  1199           _thread[u] = u + 1;

  1200           _rev_thread[u + 1] = u;

  1201           _succ_num[u] = 1;

  1202           _last_succ[u] = u;

  1203           if (_supply[u] <= 0) {

  1204             _pred_dir[u] = DIR_DOWN;

  1205             _pi[u] = 0;

  1206             _pred[u] = e;

  1207             _source[e] = _root;

  1208             _target[e] = u;

  1209             _cap[e] = INF;

  1210             _flow[e] = -_supply[u];

  1211             _cost[e] = 0;

  1212             _state[e] = STATE_TREE;

  1213           } else {

  1214             _pred_dir[u] = DIR_UP;

  1215             _pi[u] = -ART_COST;

  1216             _pred[u] = f;

  1217             _source[f] = u;

  1218             _target[f] = _root;

  1219             _cap[f] = INF;

  1220             _flow[f] = _supply[u];

  1221             _state[f] = STATE_TREE;

  1222             _cost[f] = ART_COST;

  1223             _source[e] = _root;

  1224             _target[e] = u;

  1225             _cap[e] = INF;

  1226             _flow[e] = 0;

  1227             _cost[e] = 0;

  1228             _state[e] = STATE_LOWER;

  1229             ++f;

  1230           }

  1231         }

  1232         _all_arc_num = f;

  1233       }

  1234

  1235       return true;

  1236     }

  1237

  1238     // Check if the upper bound is greater than or equal to the lower bound

  1239     // on each arc.

  1240     bool checkBoundMaps() {

  1241       for (int j = 0; j != _arc_num; ++j) {

  1242         if (_upper[j] < _lower[j]) return false;

  1243       }

  1244       return true;

  1245     }

  1246

  1247     // Find the join node

  1248     void findJoinNode() {

  1249       int u = _source[in_arc];

  1250       int v = _target[in_arc];

  1251       while (u != v) {

  1252         if (_succ_num[u] < _succ_num[v]) {

  1253           u = _parent[u];

  1254         } else {

  1255           v = _parent[v];

  1256         }

  1257       }

  1258       join = u;

  1259     }

  1260

  1261     // Find the leaving arc of the cycle and returns true if the

  1262     // leaving arc is not the same as the entering arc

  1263     bool findLeavingArc() {

  1264       // Initialize first and second nodes according to the direction

  1265       // of the cycle

  1266       int first, second;

  1267       if (_state[in_arc] == STATE_LOWER) {

  1268         first  = _source[in_arc];

  1269         second = _target[in_arc];

  1270       } else {

  1271         first  = _target[in_arc];

  1272         second = _source[in_arc];

  1273       }

  1274       delta = _cap[in_arc];

  1275       int result = 0;

  1276       Value c, d;

  1277       int e;

  1278

  1279       // Search the cycle form the first node to the join node

  1280       for (int u = first; u != join; u = _parent[u]) {

  1281         e = _pred[u];

  1282         d = _flow[e];

  1283         if (_pred_dir[u] == DIR_DOWN) {

  1284           c = _cap[e];

  1285           d = c >= MAX ? INF : c - d;

  1286         }

  1287         if (d < delta) {

  1288           delta = d;

  1289           u_out = u;

  1290           result = 1;

  1291         }

  1292       }

  1293

  1294       // Search the cycle form the second node to the join node

  1295       for (int u = second; u != join; u = _parent[u]) {

  1296         e = _pred[u];

  1297         d = _flow[e];

  1298         if (_pred_dir[u] == DIR_UP) {

  1299           c = _cap[e];

  1300           d = c >= MAX ? INF : c - d;

  1301         }

  1302         if (d <= delta) {

  1303           delta = d;

  1304           u_out = u;

  1305           result = 2;

  1306         }

  1307       }

  1308

  1309       if (result == 1) {

  1310         u_in = first;

  1311         v_in = second;

  1312       } else {

  1313         u_in = second;

  1314         v_in = first;

  1315       }

  1316       return result != 0;

  1317     }

  1318

  1319     // Change _flow and _state vectors

  1320     void changeFlow(bool change) {

  1321       // Augment along the cycle

  1322       if (delta > 0) {

  1323         Value val = _state[in_arc] * delta;

  1324         _flow[in_arc] += val;

  1325         for (int u = _source[in_arc]; u != join; u = _parent[u]) {

  1326           _flow[_pred[u]] -= _pred_dir[u] * val;

  1327         }

  1328         for (int u = _target[in_arc]; u != join; u = _parent[u]) {

  1329           _flow[_pred[u]] += _pred_dir[u] * val;

  1330         }

  1331       }

  1332       // Update the state of the entering and leaving arcs

  1333       if (change) {

  1334         _state[in_arc] = STATE_TREE;

  1335         _state[_pred[u_out]] =

  1336           (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;

  1337       } else {

  1338         _state[in_arc] = -_state[in_arc];

  1339       }

  1340     }

  1341

  1342     // Update the tree structure

  1343     void updateTreeStructure() {

  1344       int old_rev_thread = _rev_thread[u_out];

  1345       int old_succ_num = _succ_num[u_out];

  1346       int old_last_succ = _last_succ[u_out];

  1347       v_out = _parent[u_out];

  1348

  1349       // Check if u_in and u_out coincide

  1350       if (u_in == u_out) {

  1351         // Update _parent, _pred, _pred_dir

  1352         _parent[u_in] = v_in;

  1353         _pred[u_in] = in_arc;

  1354         _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;

  1355

  1356         // Update _thread and _rev_thread

  1357         if (_thread[v_in] != u_out) {

  1358           int after = _thread[old_last_succ];

  1359           _thread[old_rev_thread] = after;

  1360           _rev_thread[after] = old_rev_thread;

  1361           after = _thread[v_in];

  1362           _thread[v_in] = u_out;

  1363           _rev_thread[u_out] = v_in;

  1364           _thread[old_last_succ] = after;

  1365           _rev_thread[after] = old_last_succ;

  1366         }

  1367       } else {

  1368         // Handle the case when old_rev_thread equals to v_in

  1369         // (it also means that join and v_out coincide)

  1370         int thread_continue = old_rev_thread == v_in ?

  1371           _thread[old_last_succ] : _thread[v_in];

  1372

  1373         // Update _thread and _parent along the stem nodes (i.e. the nodes

  1374         // between u_in and u_out, whose parent have to be changed)

  1375         int stem = u_in;              // the current stem node

  1376         int par_stem = v_in;          // the new parent of stem

  1377         int next_stem;                // the next stem node

  1378         int last = _last_succ[u_in];  // the last successor of stem

  1379         int before, after = _thread[last];

  1380         _thread[v_in] = u_in;

  1381         _dirty_revs.clear();

  1382         _dirty_revs.push_back(v_in);

  1383         while (stem != u_out) {

  1384           // Insert the next stem node into the thread list

  1385           next_stem = _parent[stem];

  1386           _thread[last] = next_stem;

  1387           _dirty_revs.push_back(last);

  1388

  1389           // Remove the subtree of stem from the thread list

  1390           before = _rev_thread[stem];

  1391           _thread[before] = after;

  1392           _rev_thread[after] = before;

  1393

  1394           // Change the parent node and shift stem nodes

  1395           _parent[stem] = par_stem;

  1396           par_stem = stem;

  1397           stem = next_stem;

  1398

  1399           // Update last and after

  1400           last = _last_succ[stem] == _last_succ[par_stem] ?

  1401             _rev_thread[par_stem] : _last_succ[stem];

  1402           after = _thread[last];

  1403         }

  1404         _parent[u_out] = par_stem;

  1405         _thread[last] = thread_continue;

  1406         _rev_thread[thread_continue] = last;

  1407         _last_succ[u_out] = last;

  1408

  1409         // Remove the subtree of u_out from the thread list except for

  1410         // the case when old_rev_thread equals to v_in

  1411         if (old_rev_thread != v_in) {

  1412           _thread[old_rev_thread] = after;

  1413           _rev_thread[after] = old_rev_thread;

  1414         }

  1415

  1416         // Update _rev_thread using the new _thread values

  1417         for (int i = 0; i != int(_dirty_revs.size()); ++i) {

  1418           int u = _dirty_revs[i];

  1419           _rev_thread[_thread[u]] = u;

  1420         }

  1421

  1422         // Update _pred, _pred_dir, _last_succ and _succ_num for the

  1423         // stem nodes from u_out to u_in

  1424         int tmp_sc = 0, tmp_ls = _last_succ[u_out];

  1425         for (int u = u_out, p = _parent[u]; u != u_in; u = p, p = _parent[u]) {

  1426           _pred[u] = _pred[p];

  1427           _pred_dir[u] = -_pred_dir[p];

  1428           tmp_sc += _succ_num[u] - _succ_num[p];

  1429           _succ_num[u] = tmp_sc;

  1430           _last_succ[p] = tmp_ls;

  1431         }

  1432         _pred[u_in] = in_arc;

  1433         _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;

  1434         _succ_num[u_in] = old_succ_num;

  1435       }

  1436

  1437       // Update _last_succ from v_in towards the root

  1438       int up_limit_out = _last_succ[join] == v_in ? join : -1;

  1439       int last_succ_out = _last_succ[u_out];

  1440       for (int u = v_in; u != -1 && _last_succ[u] == v_in; u = _parent[u]) {

  1441         _last_succ[u] = last_succ_out;

  1442       }

  1443

  1444       // Update _last_succ from v_out towards the root

  1445       if (join != old_rev_thread && v_in != old_rev_thread) {

  1446         for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;

  1447              u = _parent[u]) {

  1448           _last_succ[u] = old_rev_thread;

  1449         }

  1450       }

  1451       else if (last_succ_out != old_last_succ) {

  1452         for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;

  1453              u = _parent[u]) {

  1454           _last_succ[u] = last_succ_out;

  1455         }

  1456       }

  1457

  1458       // Update _succ_num from v_in to join

  1459       for (int u = v_in; u != join; u = _parent[u]) {

  1460         _succ_num[u] += old_succ_num;

  1461       }

  1462       // Update _succ_num from v_out to join

  1463       for (int u = v_out; u != join; u = _parent[u]) {

  1464         _succ_num[u] -= old_succ_num;

  1465       }

  1466     }

  1467

  1468     // Update potentials in the subtree that has been moved

  1469     void updatePotential() {

  1470       Cost sigma = _pi[v_in] - _pi[u_in] -

  1471                    _pred_dir[u_in] * _cost[in_arc];

  1472       int end = _thread[_last_succ[u_in]];

  1473       for (int u = u_in; u != end; u = _thread[u]) {

  1474         _pi[u] += sigma;

  1475       }

  1476     }

  1477

  1478     // Heuristic initial pivots

  1479     bool initialPivots() {

  1480       Value curr, total = 0;

  1481       std::vector<Node> supply_nodes, demand_nodes;

  1482       for (NodeIt u(_graph); u != INVALID; ++u) {

  1483         curr = _supply[_node_id[u]];

  1484         if (curr > 0) {

  1485           total += curr;

  1486           supply_nodes.push_back(u);

  1487         }

  1488         else if (curr < 0) {

  1489           demand_nodes.push_back(u);

  1490         }

  1491       }

  1492       if (_sum_supply > 0) total -= _sum_supply;

  1493       if (total <= 0) return true;

  1494

  1495       IntVector arc_vector;

  1496       if (_sum_supply >= 0) {

  1497         if (supply_nodes.size() == 1 && demand_nodes.size() == 1) {

  1498           // Perform a reverse graph search from the sink to the source

  1499           typename GR::template NodeMap<bool> reached(_graph, false);

  1500           Node s = supply_nodes[0], t = demand_nodes[0];

  1501           std::vector<Node> stack;

  1502           reached[t] = true;

  1503           stack.push_back(t);

  1504           while (!stack.empty()) {

  1505             Node u, v = stack.back();

  1506             stack.pop_back();

  1507             if (v == s) break;

  1508             for (InArcIt a(_graph, v); a != INVALID; ++a) {

  1509               if (reached[u = _graph.source(a)]) continue;

  1510               int j = _arc_id[a];

  1511               if (_cap[j] >= total) {

  1512                 arc_vector.push_back(j);

  1513                 reached[u] = true;

  1514                 stack.push_back(u);

  1515               }

  1516             }

  1517           }

  1518         } else {

  1519           // Find the min. cost incoming arc for each demand node

  1520           for (int i = 0; i != int(demand_nodes.size()); ++i) {

  1521             Node v = demand_nodes[i];

  1522             Cost c, min_cost = std::numeric_limits<Cost>::max();

  1523             Arc min_arc = INVALID;

  1524             for (InArcIt a(_graph, v); a != INVALID; ++a) {

  1525               c = _cost[_arc_id[a]];

  1526               if (c < min_cost) {

  1527                 min_cost = c;

  1528                 min_arc = a;

  1529               }

  1530             }

  1531             if (min_arc != INVALID) {

  1532               arc_vector.push_back(_arc_id[min_arc]);

  1533             }

  1534           }

  1535         }

  1536       } else {

  1537         // Find the min. cost outgoing arc for each supply node

  1538         for (int i = 0; i != int(supply_nodes.size()); ++i) {

  1539           Node u = supply_nodes[i];

  1540           Cost c, min_cost = std::numeric_limits<Cost>::max();

  1541           Arc min_arc = INVALID;

  1542           for (OutArcIt a(_graph, u); a != INVALID; ++a) {

  1543             c = _cost[_arc_id[a]];

  1544             if (c < min_cost) {

  1545               min_cost = c;

  1546               min_arc = a;

  1547             }

  1548           }

  1549           if (min_arc != INVALID) {

  1550             arc_vector.push_back(_arc_id[min_arc]);

  1551           }

  1552         }

  1553       }

  1554

  1555       // Perform heuristic initial pivots

  1556       for (int i = 0; i != int(arc_vector.size()); ++i) {

  1557         in_arc = arc_vector[i];

  1558         if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] -

  1559             _pi[_target[in_arc]]) >= 0) continue;

  1560         findJoinNode();

  1561         bool change = findLeavingArc();

  1562         if (delta >= MAX) return false;

  1563         changeFlow(change);

  1564         if (change) {

  1565           updateTreeStructure();

  1566           updatePotential();

  1567         }

  1568       }

  1569       return true;

  1570     }

  1571

  1572     // Execute the algorithm

  1573     ProblemType start(PivotRule pivot_rule) {

  1574       // Select the pivot rule implementation

  1575       switch (pivot_rule) {

  1576         case FIRST_ELIGIBLE:

  1577           return start<FirstEligiblePivotRule>();

  1578         case BEST_ELIGIBLE:

  1579           return start<BestEligiblePivotRule>();

  1580         case BLOCK_SEARCH:

  1581           return start<BlockSearchPivotRule>();

  1582         case CANDIDATE_LIST:

  1583           return start<CandidateListPivotRule>();

  1584         case ALTERING_LIST:

  1585           return start<AlteringListPivotRule>();

  1586       }

  1587       return INFEASIBLE; // avoid warning

  1588     }

  1589

  1590     template <typename PivotRuleImpl>

  1591     ProblemType start() {

  1592       PivotRuleImpl pivot(*this);

  1593

  1594       // Perform heuristic initial pivots

  1595       if (!initialPivots()) return UNBOUNDED;

  1596

  1597       // Execute the Network Simplex algorithm

  1598       while (pivot.findEnteringArc()) {

  1599         findJoinNode();

  1600         bool change = findLeavingArc();

  1601         if (delta >= MAX) return UNBOUNDED;

  1602         changeFlow(change);

  1603         if (change) {

  1604           updateTreeStructure();

  1605           updatePotential();

  1606         }

  1607       }

  1608

  1609       // Check feasibility

  1610       for (int e = _search_arc_num; e != _all_arc_num; ++e) {

  1611         if (_flow[e] != 0) return INFEASIBLE;

  1612       }

  1613

  1614       // Transform the solution and the supply map to the original form

  1615       if (_has_lower) {

  1616         for (int i = 0; i != _arc_num; ++i) {

  1617           Value c = _lower[i];

  1618           if (c != 0) {

  1619             _flow[i] += c;

  1620             _supply[_source[i]] += c;

  1621             _supply[_target[i]] -= c;

  1622           }

  1623         }

  1624       }

  1625

  1626       // Shift potentials to meet the requirements of the GEQ/LEQ type

  1627       // optimality conditions

  1628       if (_sum_supply == 0) {

  1629         if (_stype == GEQ) {

  1630           Cost max_pot = -std::numeric_limits<Cost>::max();

  1631           for (int i = 0; i != _node_num; ++i) {

  1632             if (_pi[i] > max_pot) max_pot = _pi[i];

  1633           }

  1634           if (max_pot > 0) {

  1635             for (int i = 0; i != _node_num; ++i)

  1636               _pi[i] -= max_pot;

  1637           }

  1638         } else {

  1639           Cost min_pot = std::numeric_limits<Cost>::max();

  1640           for (int i = 0; i != _node_num; ++i) {

  1641             if (_pi[i] < min_pot) min_pot = _pi[i];

  1642           }

  1643           if (min_pot < 0) {

  1644             for (int i = 0; i != _node_num; ++i)

  1645               _pi[i] -= min_pot;

  1646           }

  1647         }

  1648       }

  1649

  1650       return OPTIMAL;

  1651     }

  1652

  1653   }; //class NetworkSimplex

  1654

  1655   ///@}

  1656

  1657 } //namespace lemon

  1658

  1659 #endif //LEMON_NETWORK_SIMPLEX_H