#define LEMON_NETWORK_SIMPLEX_H

 #define LEMON_NETWORK_SIMPLEX_H

 /// \ingroup min_cost_flow

 /// \ingroup min_cost_flow

 ///

 ///

 /// \file

 /// \file

 #include <vector>

 #include <vector>

 #include <limits>

 #include <limits>

 #include <algorithm>

 #include <algorithm>

 namespace lemon {

 namespace lemon {

   /// \addtogroup min_cost_flow

   /// \addtogroup min_cost_flow

   /// @{

   /// @{

   /// for finding a \ref min_cost_flow "minimum cost flow".

   /// for finding a \ref min_cost_flow "minimum cost flow".

   ///

   ///

   /// for finding a \ref min_cost_flow "minimum cost flow".

   /// for finding a \ref min_cost_flow "minimum cost flow".

   ///

   ///

   ///

   ///

   ///

   ///

   class NetworkSimplex

   class NetworkSimplex

   {

   {

     typedef std::vector<Arc> ArcVector;

     typedef std::vector<Arc> ArcVector;

     typedef std::vector<Node> NodeVector;

     typedef std::vector<Node> NodeVector;

     typedef std::vector<int> IntVector;

     typedef std::vector<int> IntVector;

     typedef std::vector<bool> BoolVector;

     typedef std::vector<bool> BoolVector;

     // State constants for arcs

     // State constants for arcs

     enum ArcStateEnum {

     enum ArcStateEnum {

       STATE_UPPER = -1,

       STATE_UPPER = -1,

       STATE_TREE  =  0,

       STATE_TREE  =  0,

       STATE_LOWER =  1

       STATE_LOWER =  1

     };

     };

   private:

   private:

     // Result maps

     // Result maps

     FlowMap *_flow_map;

     FlowMap *_flow_map;

     PotentialMap *_potential_map;

     PotentialMap *_potential_map;

     bool _local_flow;

     bool _local_flow;

     bool _local_potential;

     bool _local_potential;

     IntNodeMap _node_id;

     IntNodeMap _node_id;

     ArcVector _arc_ref;

     ArcVector _arc_ref;

     IntVector _source;

     IntVector _source;

     IntVector _target;

     IntVector _target;

     // Data for storing the spanning tree structure

     // Data for storing the spanning tree structure

     IntVector _parent;

     IntVector _parent;

     IntVector _pred;

     IntVector _pred;

     IntVector _thread;

     IntVector _thread;

     // Temporary data used in the current pivot iteration

     // Temporary data used in the current pivot iteration

     int in_arc, join, u_in, v_in, u_out, v_out;

     int in_arc, join, u_in, v_in, u_out, v_out;

     int first, second, right, last;

     int first, second, right, last;

     int stem, par_stem, new_stem;

     int stem, par_stem, new_stem;

   private:

   private:

     class FirstEligiblePivotRule

     class FirstEligiblePivotRule

     {

     {

     private:

     private:

       // References to the NetworkSimplex class

       // References to the NetworkSimplex class

       const IntVector  &_source;

       const IntVector  &_source;

       const IntVector  &_target;

       const IntVector  &_target;

       const IntVector  &_state;

       const IntVector  &_state;

       int &_in_arc;

       int &_in_arc;

       int _arc_num;

       int _arc_num;

       // Pivot rule data

       // Pivot rule data

       int _next_arc;

       int _next_arc;

     public:

     public:

       FirstEligiblePivotRule(NetworkSimplex &ns) :

       FirstEligiblePivotRule(NetworkSimplex &ns) :

         _source(ns._source), _target(ns._target),

         _source(ns._source), _target(ns._target),

         _cost(ns._cost), _state(ns._state), _pi(ns._pi),

         _cost(ns._cost), _state(ns._state), _pi(ns._pi),

         _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)

         _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)

       {}

       {}

       bool findEnteringArc() {

       bool findEnteringArc() {

         for (int e = _next_arc; e < _arc_num; ++e) {

         for (int e = _next_arc; e < _arc_num; ++e) {

           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

           if (c < 0) {

           if (c < 0) {

             _in_arc = e;

             _in_arc = e;

             _next_arc = e + 1;

             _next_arc = e + 1;

       }

       }

     }; //class FirstEligiblePivotRule

     }; //class FirstEligiblePivotRule

     class BestEligiblePivotRule

     class BestEligiblePivotRule

     {

     {

     private:

     private:

       // References to the NetworkSimplex class

       // References to the NetworkSimplex class

       const IntVector  &_source;

       const IntVector  &_source;

       const IntVector  &_target;

       const IntVector  &_target;

       const IntVector  &_state;

       const IntVector  &_state;

       int &_in_arc;

       int &_in_arc;

       int _arc_num;

       int _arc_num;

     public:

     public:

       BestEligiblePivotRule(NetworkSimplex &ns) :

       BestEligiblePivotRule(NetworkSimplex &ns) :

         _source(ns._source), _target(ns._target),

         _source(ns._source), _target(ns._target),

         _cost(ns._cost), _state(ns._state), _pi(ns._pi),

         _cost(ns._cost), _state(ns._state), _pi(ns._pi),

         _in_arc(ns.in_arc), _arc_num(ns._arc_num)

         _in_arc(ns.in_arc), _arc_num(ns._arc_num)

       {}

       {}

       bool findEnteringArc() {

       bool findEnteringArc() {

         for (int e = 0; e < _arc_num; ++e) {

         for (int e = 0; e < _arc_num; ++e) {

           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

           if (c < min) {

           if (c < min) {

             min = c;

             min = c;

             _in_arc = e;

             _in_arc = e;

       }

       }

     }; //class BestEligiblePivotRule

     }; //class BestEligiblePivotRule

     class BlockSearchPivotRule

     class BlockSearchPivotRule

     {

     {

     private:

     private:

       // References to the NetworkSimplex class

       // References to the NetworkSimplex class

       const IntVector  &_source;

       const IntVector  &_source;

       const IntVector  &_target;

       const IntVector  &_target;

       const IntVector  &_state;

       const IntVector  &_state;

       int &_in_arc;

       int &_in_arc;

       int _arc_num;

       int _arc_num;

       // Pivot rule data

       // Pivot rule data

       int _block_size;

       int _block_size;

       int _next_arc;

       int _next_arc;

     public:

     public:

       BlockSearchPivotRule(NetworkSimplex &ns) :

       BlockSearchPivotRule(NetworkSimplex &ns) :

         _source(ns._source), _target(ns._target),

         _source(ns._source), _target(ns._target),

         _cost(ns._cost), _state(ns._state), _pi(ns._pi),

         _cost(ns._cost), _state(ns._state), _pi(ns._pi),

         _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)

         _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)

       {

       {

         _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),

         _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),

                                 MIN_BLOCK_SIZE );

                                 MIN_BLOCK_SIZE );

       }

       }

       bool findEnteringArc() {

       bool findEnteringArc() {

         int cnt = _block_size;

         int cnt = _block_size;

         int e, min_arc = _next_arc;

         int e, min_arc = _next_arc;

         for (e = _next_arc; e < _arc_num; ++e) {

         for (e = _next_arc; e < _arc_num; ++e) {

           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

           if (c < min) {

           if (c < min) {

       }

       }

     }; //class BlockSearchPivotRule

     }; //class BlockSearchPivotRule

     class CandidateListPivotRule

     class CandidateListPivotRule

     {

     {

     private:

     private:

       // References to the NetworkSimplex class

       // References to the NetworkSimplex class

       const IntVector  &_source;

       const IntVector  &_source;

       const IntVector  &_target;

       const IntVector  &_target;

       const IntVector  &_state;

       const IntVector  &_state;

       int &_in_arc;

       int &_in_arc;

       int _arc_num;

       int _arc_num;

       // Pivot rule data

       // Pivot rule data

       IntVector _candidates;

       IntVector _candidates;

         _candidates.resize(_list_length);

         _candidates.resize(_list_length);

       }

       }

       /// Find next entering arc

       /// Find next entering arc

       bool findEnteringArc() {

       bool findEnteringArc() {

         int e, min_arc = _next_arc;

         int e, min_arc = _next_arc;

         if (_curr_length > 0 && _minor_count < _minor_limit) {

         if (_curr_length > 0 && _minor_count < _minor_limit) {

           // Minor iteration: select the best eligible arc from the

           // Minor iteration: select the best eligible arc from the

           // current candidate list

           // current candidate list

           ++_minor_count;

           ++_minor_count;

       }

       }

     }; //class CandidateListPivotRule

     }; //class CandidateListPivotRule

     class AlteringListPivotRule

     class AlteringListPivotRule

     {

     {

     private:

     private:

       // References to the NetworkSimplex class

       // References to the NetworkSimplex class

       const IntVector  &_source;

       const IntVector  &_source;

       const IntVector  &_target;

       const IntVector  &_target;

       const IntVector  &_state;

       const IntVector  &_state;

       int &_in_arc;

       int &_in_arc;

       int _arc_num;

       int _arc_num;

       // Pivot rule data

       // Pivot rule data

       int _block_size, _head_length, _curr_length;

       int _block_size, _head_length, _curr_length;

       int _next_arc;

       int _next_arc;

       IntVector _candidates;

       IntVector _candidates;

       // Functor class to compare arcs during sort of the candidate list

       // Functor class to compare arcs during sort of the candidate list

       class SortFunc

       class SortFunc

       {

       {

       private:

       private:

       public:

       public:

         bool operator()(int left, int right) {

         bool operator()(int left, int right) {

           return _map[left] > _map[right];

           return _map[left] > _map[right];

         }

         }

       };

       };

       SortFunc _sort_func;

       SortFunc _sort_func;

     public:

     public:

       AlteringListPivotRule(NetworkSimplex &ns) :

       AlteringListPivotRule(NetworkSimplex &ns) :

         _source(ns._source), _target(ns._target),

         _source(ns._source), _target(ns._target),

         _cost(ns._cost), _state(ns._state), _pi(ns._pi),

         _cost(ns._cost), _state(ns._state), _pi(ns._pi),

         _in_arc(ns.in_arc), _arc_num(ns._arc_num),

         _in_arc(ns.in_arc), _arc_num(ns._arc_num),

         _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)

         _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)

                                  MIN_HEAD_LENGTH );

                                  MIN_HEAD_LENGTH );

         _candidates.resize(_head_length + _block_size);

         _candidates.resize(_head_length + _block_size);

         _curr_length = 0;

         _curr_length = 0;

       }

       }

       bool findEnteringArc() {

       bool findEnteringArc() {

         // Check the current candidate list

         // Check the current candidate list

         int e;

         int e;

         for (int i = 0; i < _curr_length; ++i) {

         for (int i = 0; i < _curr_length; ++i) {

           e = _candidates[i];

           e = _candidates[i];

     }; //class AlteringListPivotRule

     }; //class AlteringListPivotRule

   public:

   public:

     ///

     ///

     /// \param graph The digraph the algorithm runs on.

     /// \param graph The digraph the algorithm runs on.

       _flow_map(NULL), _potential_map(NULL),

       _flow_map(NULL), _potential_map(NULL),

       _local_flow(false), _local_potential(false),

       _local_flow(false), _local_potential(false),

       _node_id(graph)

       _node_id(graph)

     /// Destructor.

     /// Destructor.

     ~NetworkSimplex() {

     ~NetworkSimplex() {

       if (_local_flow) delete _flow_map;

       if (_local_flow) delete _flow_map;

       if (_local_potential) delete _potential_map;

       if (_local_potential) delete _potential_map;

     }

     }

     /// \brief Set the flow map.

     /// \brief Set the flow map.

     ///

     ///

     /// This function sets the flow map.

     /// This function sets the flow map.

     ///

     ///

     /// \return <tt>(*this)</tt>

     /// \return <tt>(*this)</tt>

       if (_local_flow) {

       if (_local_flow) {

         delete _flow_map;

         delete _flow_map;

         _local_flow = false;

         _local_flow = false;

       }

       }

       _flow_map = ↦

       _flow_map = ↦

       return *this;

       return *this;

     }

     }

     /// \brief Set the potential map.

     /// \brief Set the potential map.

     ///

     ///

     ///

     ///

     /// \return <tt>(*this)</tt>

     /// \return <tt>(*this)</tt>

       if (_local_potential) {

       if (_local_potential) {

         delete _potential_map;

         delete _potential_map;

         _local_potential = false;

         _local_potential = false;

       }

       }

       _potential_map = ↦

       _potential_map = ↦

       return *this;

       return *this;

     }

     }

     /// @{

     /// @{

     /// \brief Run the algorithm.

     /// \brief Run the algorithm.

     ///

     ///

     /// This function runs the algorithm.

     /// This function runs the algorithm.

     ///

     ///

     /// \return \c true if a feasible flow can be found.

     /// \return \c true if a feasible flow can be found.

       return init() && start(pivot_rule);

       return init() && start(pivot_rule);

     }

     }

     /// @}

     /// @}

     /// \name Query Functions

     /// \name Query Functions

     /// The results of the algorithm can be obtained using these

     /// The results of the algorithm can be obtained using these

     /// functions.\n

     /// functions.\n

     /// @{

     /// @{

     /// \brief Return a const reference to the flow map.

     /// \brief Return a const reference to the flow map.

     ///

     ///

     /// This function returns a const reference to an arc map storing

     /// This function returns a const reference to an arc map storing

     /// the found flow.

     /// the found flow.

     /// \pre \ref run() must be called before using this function.

     /// \pre \ref run() must be called before using this function.

     const FlowMap& flowMap() const {

     const FlowMap& flowMap() const {

       return *_flow_map;

       return *_flow_map;

     }

     }

     /// \brief Return a const reference to the potential map

     /// \brief Return a const reference to the potential map

     /// (the dual solution).

     /// (the dual solution).

     ///

     ///

     /// This function returns a const reference to a node map storing

     /// This function returns a const reference to a node map storing

     ///

     ///

     /// \pre \ref run() must be called before using this function.

     /// \pre \ref run() must be called before using this function.

     const PotentialMap& potentialMap() const {

     const PotentialMap& potentialMap() const {

       return *_potential_map;

       return *_potential_map;

     }

     }

     /// @}

     /// @}

   private:

   private:

       // Initialize vectors

       // Initialize vectors

       _node_num = countNodes(_graph);

       _node_num = countNodes(_graph);

       _arc_num = countArcs(_graph);

       _arc_num = countArcs(_graph);

       int all_node_num = _node_num + 1;

       int all_node_num = _node_num + 1;

       int all_arc_num = _arc_num + _node_num;

       int all_arc_num = _arc_num + _node_num;

       _arc_ref.resize(_arc_num);

       _arc_ref.resize(_arc_num);

       _source.resize(all_arc_num);

       _source.resize(all_arc_num);

       _target.resize(all_arc_num);

       _target.resize(all_arc_num);

       _last_succ.resize(all_node_num);

       _last_succ.resize(all_node_num);

       _state.resize(all_arc_num, STATE_LOWER);

       _state.resize(all_arc_num, STATE_LOWER);

       // Initialize node related data

       // Initialize node related data

       bool valid_supply = true;

       bool valid_supply = true;

         int i = 0;

         int i = 0;

         for (NodeIt n(_graph); n != INVALID; ++n, ++i) {

         for (NodeIt n(_graph); n != INVALID; ++n, ++i) {

           _node_id[n] = i;

           _node_id[n] = i;

           sum += _supply[i];

           sum += _supply[i];

         }

         }

         valid_supply = (sum == 0);

         valid_supply = (sum == 0);

       } else {

       } else {

         int i = 0;

         int i = 0;

         for (NodeIt n(_graph); n != INVALID; ++n, ++i) {

         for (NodeIt n(_graph); n != INVALID; ++n, ++i) {

           _node_id[n] = i;

           _node_id[n] = i;

           _supply[i] = 0;

           _supply[i] = 0;

         }

         }

       }

       }

       if (!valid_supply) return false;

       if (!valid_supply) return false;

       // Set data for the artificial root node

       // Set data for the artificial root node

       _root = _node_num;

       _root = _node_num;

         _arc_ref[i] = e;

         _arc_ref[i] = e;

         if ((i += k) >= _arc_num) i = (i % k) + 1;

         if ((i += k) >= _arc_num) i = (i % k) + 1;

       }

       }

       // Initialize arc maps

       // Initialize arc maps

       }

       }

       // Remove non-zero lower bounds

       // Remove non-zero lower bounds

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

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

           if (c != 0) {

           if (c != 0) {

             _cap[i] -= c;

             _cap[i] -= c;

             _supply[_source[i]] -= c;

             _supply[_source[i]] -= c;

             _supply[_target[i]] += c;

             _supply[_target[i]] += c;

           }

           }

         }

         }

       }

       }

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

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

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

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

         _thread[u] = u + 1;

         _thread[u] = u + 1;

         _rev_thread[u + 1] = u;

         _rev_thread[u + 1] = u;

         _succ_num[u] = 1;

         _succ_num[u] = 1;

         _last_succ[u] = u;

         _last_succ[u] = u;

         first  = _target[in_arc];

         first  = _target[in_arc];

         second = _source[in_arc];

         second = _source[in_arc];

       }

       }

       delta = _cap[in_arc];

       delta = _cap[in_arc];

       int result = 0;

       int result = 0;

       int e;

       int e;

       // Search the cycle along the path form the first node to the root

       // Search the cycle along the path form the first node to the root

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

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

         e = _pred[u];

         e = _pred[u];

     // Change _flow and _state vectors

     // Change _flow and _state vectors

     void changeFlow(bool change) {

     void changeFlow(bool change) {

       // Augment along the cycle

       // Augment along the cycle

       if (delta > 0) {

       if (delta > 0) {

         _flow[in_arc] += val;

         _flow[in_arc] += val;

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

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

           _flow[_pred[u]] += _forward[u] ? -val : val;

           _flow[_pred[u]] += _forward[u] ? -val : val;

         }

         }

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

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

       }

       }

     }

     }

     // Update potentials

     // Update potentials

     void updatePotential() {

     void updatePotential() {

         _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :

         _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :

         _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];

         _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];

       if (_succ_num[u_in] > _node_num / 2) {

       if (_succ_num[u_in] > _node_num / 2) {

         // Update in the upper subtree (which contains the root)

         // Update in the upper subtree (which contains the root)

         int before = _rev_thread[u_in];

         int before = _rev_thread[u_in];

         }

         }

       }

       }

     }

     }

     // Execute the algorithm

     // Execute the algorithm

       // Select the pivot rule implementation

       // Select the pivot rule implementation

       switch (pivot_rule) {

       switch (pivot_rule) {

           return start<FirstEligiblePivotRule>();

           return start<FirstEligiblePivotRule>();

           return start<BestEligiblePivotRule>();

           return start<BestEligiblePivotRule>();

           return start<BlockSearchPivotRule>();

           return start<BlockSearchPivotRule>();

           return start<CandidateListPivotRule>();

           return start<CandidateListPivotRule>();

           return start<AlteringListPivotRule>();

           return start<AlteringListPivotRule>();

       }

       }

       return false;

       return false;

     }

     }

     bool start() {

     bool start() {

       while (pivot.findEnteringArc()) {

       while (pivot.findEnteringArc()) {

         findJoinNode();

         findJoinNode();

         bool change = findLeavingArc();

         bool change = findLeavingArc();

         changeFlow(change);

         changeFlow(change);

         if (change) {

         if (change) {

       for (int e = _arc_num; e != _arc_num + _node_num; ++e) {

       for (int e = _arc_num; e != _arc_num + _node_num; ++e) {

         if (_flow[e] > 0) return false;

         if (_flow[e] > 0) return false;

       }

       }

       // Copy flow values to _flow_map

       // Copy flow values to _flow_map

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

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

           Arc e = _arc_ref[i];

           Arc e = _arc_ref[i];

         }

         }

       } else {

       } else {

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

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

           _flow_map->set(_arc_ref[i], _flow[i]);

           _flow_map->set(_arc_ref[i], _flow[i]);

         }

         }