/// \brief The network simplex algorithm for finding a minimum cost flow.

 /// \brief Network simplex algorithm for finding a minimum cost flow.

 #include <vector>

 #include <lemon/math.h>

 /// \brief The pivot rule used in the algorithm.

 //#define FIRST_ELIGIBLE_PIVOT

 //#define BEST_ELIGIBLE_PIVOT

 #define EDGE_BLOCK_PIVOT

 //#define CANDIDATE_LIST_PIVOT

 //#define SORTED_LIST_PIVOT

 //#define _DEBUG_ITER_

 // State constant for edges at their lower bounds.

 #define LOWER   1

 // State constant for edges in the spanning tree.

 #define TREE    0

 // State constant for edges at their upper bounds.

 #define UPPER   -1

 #ifdef EDGE_BLOCK_PIVOT

   #include <lemon/math.h>

   #define MIN_BLOCK_SIZE        10

 #endif

 #ifdef CANDIDATE_LIST_PIVOT

   #include <vector>

   #define LIST_LENGTH_DIV       50

   #define MINOR_LIMIT_DIV       200

 #endif

 #ifdef SORTED_LIST_PIVOT

   #include <vector>

   #include <algorithm>

   #define LIST_LENGTH_DIV       100

   #define LOWER_DIV             4

 #endif

   /// \param Graph The directed graph type the algorithm runs on.

   /// \tparam Graph The directed graph type the algorithm runs on.

   /// \param LowerMap The type of the lower bound map.

   /// \tparam LowerMap The type of the lower bound map.

   /// \param CapacityMap The type of the capacity (upper bound) map.

   /// \tparam CapacityMap The type of the capacity (upper bound) map.

   /// \param CostMap The type of the cost (length) map.

   /// \tparam CostMap The type of the cost (length) map.

   /// \param SupplyMap The type of the supply map.

   /// \tparam SupplyMap The type of the supply map.

   /// - Edge capacities and costs should be non-negative integers.

   /// - Edge capacities and costs should be \e non-negative \e integers.

   ///   However \c CostMap::Value should be signed type.

   /// - Supply values should be \e signed \e integers.

   /// - Supply values should be signed integers.

   /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.

   /// - \c LowerMap::Value must be convertible to

   /// - \c CapacityMap::Value and \c SupplyMap::Value must be

   ///   \c CapacityMap::Value and \c CapacityMap::Value must be

   ///   convertible to each other.

   ///   convertible to \c SupplyMap::Value.

   /// - All value types must be convertible to \c CostMap::Value, which

   ///   must be signed type.

   ///

   /// \note \ref NetworkSimplex provides six different pivot rule

   /// implementations that significantly affect the efficiency of the

   /// algorithm.

   /// By default a combined pivot rule is used, which is the fastest

   /// implementation according to our benchmark tests.

   /// Another pivot rule can be selected using \ref run() function

   /// with the proper parameter.

              typename CapacityMap = LowerMap,

              typename CapacityMap = typename Graph::template EdgeMap<int>,

              typename SupplyMap = typename Graph::template NodeMap

              typename SupplyMap = typename Graph::template NodeMap<int> >

                                   <typename CapacityMap::Value> >

   protected:

   public:

     /// Enum type to select the pivot rule used by \ref run().

     enum PivotRuleEnum {

       FIRST_ELIGIBLE_PIVOT,

       BEST_ELIGIBLE_PIVOT,

       BLOCK_SEARCH_PIVOT,

       LIMITED_SEARCH_PIVOT,

       CANDIDATE_LIST_PIVOT,

       COMBINED_PIVOT

     };

   private:

     /// \brief Map adaptor class for handling reduced edge costs.

     ///

       const SGraph &gr;

       const SGraph &_gr;

       const SCostMap &cost_map;

       const SCostMap &_cost_map;

       const SPotentialMap &pot_map;

       const SPotentialMap &_pot_map;

       ReducedCostMap( const SGraph &_gr,

       ///\e

                       const SCostMap &_cm,

       ReducedCostMap( const SGraph &gr,

                       const SPotentialMap &_pm ) :

                       const SCostMap &cost_map,

         gr(_gr), cost_map(_cm), pot_map(_pm) {}

                       const SPotentialMap &pot_map ) :

         _gr(gr), _cost_map(cost_map), _pot_map(pm) {}

       ///\e

         return cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];

         return _cost_map[e] + _pot_map[_gr.source(e)]

                            - _pot_map[_gr.target(e)];

   protected:

   private:

     /// The directed graph the algorithm runs on.

     /// \brief Implementation of the "First Eligible" pivot rule for the

     SGraph graph;

     /// \ref NetworkSimplex "network simplex" algorithm.

     /// The original graph.

     ///

     const Graph &graph_ref;

     /// This class implements the "First Eligible" pivot rule

     /// The original lower bound map.

     /// for the \ref NetworkSimplex "network simplex" algorithm.

     const LowerMap *lower;

     class FirstEligiblePivotRule

     /// The capacity map.

     {

     SCapacityMap capacity;

     private:

     /// The cost map.

     SCostMap cost;

       NetworkSimplex &_ns;

     /// The supply map.

       EdgeIt _next_edge;

     SSupplyMap supply;

     /// The reduced cost map.

     public:

     ReducedCostMap red_cost;

     bool valid_supply;

       /// Constructor.

       FirstEligiblePivotRule(NetworkSimplex &ns) :

     /// The edge map of the current flow.

         _ns(ns), _next_edge(ns._graph) {}

     SCapacityMap flow;

     /// The potential node map.

       /// Finds the next entering edge.

     SPotentialMap potential;

       bool findEnteringEdge() {

     FlowMap flow_result;

         for (EdgeIt e = _next_edge; e != INVALID; ++e) {

     PotentialMap potential_result;

           if (_ns._state[e] * _ns._red_cost[e] < 0) {

             _ns._in_edge = e;

     /// Node reference for the original graph.

             _next_edge = ++e;

     NodeRefMap node_ref;

             return true;

     /// Edge reference for the original graph.

           }

     EdgeRefMap edge_ref;

         }

         for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {

     /// The \c depth node map of the spanning tree structure.

           if (_ns._state[e] * _ns._red_cost[e] < 0) {

     IntNodeMap depth;

             _ns._in_edge = e;

     /// The \c parent node map of the spanning tree structure.

             _next_edge = ++e;

     NodeNodeMap parent;

             return true;

     /// The \c pred_edge node map of the spanning tree structure.

           }

     EdgeNodeMap pred_edge;

         }

     /// The \c thread node map of the spanning tree structure.

         return false;

     NodeNodeMap thread;

       }

     /// The \c forward node map of the spanning tree structure.

     }; //class FirstEligiblePivotRule

     BoolNodeMap forward;

     /// The state edge map.

     /// \brief Implementation of the "Best Eligible" pivot rule for the

     IntEdgeMap state;

     /// \ref NetworkSimplex "network simplex" algorithm.

     /// The root node of the starting spanning tree.

     ///

     Node root;

     /// This class implements the "Best Eligible" pivot rule

     /// for the \ref NetworkSimplex "network simplex" algorithm.

     class BestEligiblePivotRule

     {

     private:

       NetworkSimplex &_ns;

     public:

       /// Constructor.

       BestEligiblePivotRule(NetworkSimplex &ns) : _ns(ns) {}

       /// Finds the next entering edge.

       bool findEnteringEdge() {

         Cost min = 0;

         for (EdgeIt e(_ns._graph); e != INVALID; ++e) {

           if (_ns._state[e] * _ns._red_cost[e] < min) {

             min = _ns._state[e] * _ns._red_cost[e];

             _ns._in_edge = e;

           }

         }

         return min < 0;

       }

     }; //class BestEligiblePivotRule

     /// \brief Implementation of the "Block Search" pivot rule for the

     /// \ref NetworkSimplex "network simplex" algorithm.

     ///

     /// This class implements the "Block Search" pivot rule

     /// for the \ref NetworkSimplex "network simplex" algorithm.

     class BlockSearchPivotRule

     {

     private:

       NetworkSimplex &_ns;

       EdgeIt _next_edge, _min_edge;

       int _block_size;

       static const int MIN_BLOCK_SIZE = 10;

     public:

       /// Constructor.

       BlockSearchPivotRule(NetworkSimplex &ns) :

         _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)

       {

         _block_size = 2 * int(sqrt(countEdges(_ns._graph)));

         if (_block_size < MIN_BLOCK_SIZE) _block_size = MIN_BLOCK_SIZE;

       }

       /// Finds the next entering edge.

       bool findEnteringEdge() {

         Cost curr, min = 0;

         int cnt = 0;

         for (EdgeIt e = _next_edge; e != INVALID; ++e) {

           if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {

             min = curr;

             _min_edge = e;

           }

           if (++cnt == _block_size) {

             if (min < 0) break;

             cnt = 0;

           }

         }

         if (min == 0) {

           for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {

             if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {

               min = curr;

               _min_edge = e;

             }

             if (++cnt == _block_size) {

               if (min < 0) break;

               cnt = 0;

             }

           }

         }

         _ns._in_edge = _min_edge;

         _next_edge = ++_min_edge;

         return min < 0;

       }

     }; //class BlockSearchPivotRule

     /// \brief Implementation of the "Limited Search" pivot rule for the

     /// \ref NetworkSimplex "network simplex" algorithm.

     ///

     /// This class implements the "Limited Search" pivot rule

     /// for the \ref NetworkSimplex "network simplex" algorithm.

     class LimitedSearchPivotRule

     {

     private:

       NetworkSimplex &_ns;

       EdgeIt _next_edge, _min_edge;

       int _sample_size;

       static const int MIN_SAMPLE_SIZE = 10;

       static const double SAMPLE_SIZE_FACTOR = 0.0015;

     public:

       /// Constructor.

       LimitedSearchPivotRule(NetworkSimplex &ns) :

         _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)

       {

         _sample_size = int(SAMPLE_SIZE_FACTOR * countEdges(_ns._graph));

         if (_sample_size < MIN_SAMPLE_SIZE)

           _sample_size = MIN_SAMPLE_SIZE;

       }

       /// Finds the next entering edge.

       bool findEnteringEdge() {

         Cost curr, min = 0;

         int cnt = 0;

         for (EdgeIt e = _next_edge; e != INVALID; ++e) {

           if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {

             min = curr;

             _min_edge = e;

           }

           if (curr < 0 && ++cnt == _sample_size) break;

         }

         if (min == 0) {

           for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {

             if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {

               min = curr;

               _min_edge = e;

             }

             if (curr < 0 && ++cnt == _sample_size) break;

           }

         }

         _ns._in_edge = _min_edge;

         _next_edge = ++_min_edge;

         return min < 0;

       }

     }; //class LimitedSearchPivotRule

     /// \brief Implementation of the "Candidate List" pivot rule for the

     /// \ref NetworkSimplex "network simplex" algorithm.

     ///

     /// This class implements the "Candidate List" pivot rule

     /// for the \ref NetworkSimplex "network simplex" algorithm.

     class CandidateListPivotRule

     {

     private:

       NetworkSimplex &_ns;

       // The list of candidate edges.

       std::vector<Edge> _candidates;

       // The maximum length of the edge list.

       int _list_length;

       // The maximum number of minor iterations between two major

       // itarations.

       int _minor_limit;

       int _minor_count;

       EdgeIt _next_edge;

       static const double LIST_LENGTH_FACTOR = 0.002;

       static const double MINOR_LIMIT_FACTOR = 0.1;

       static const int MIN_LIST_LENGTH = 10;

       static const int MIN_MINOR_LIMIT = 2;

     public:

       /// Constructor.

       CandidateListPivotRule(NetworkSimplex &ns) :

         _ns(ns), _next_edge(ns._graph)

       {

         int edge_num = countEdges(_ns._graph);

         _minor_count = 0;

         _list_length = int(edge_num * LIST_LENGTH_FACTOR);

         if (_list_length < MIN_LIST_LENGTH)

           _list_length = MIN_LIST_LENGTH;

         _minor_limit = int(_list_length * MINOR_LIMIT_FACTOR);

         if (_minor_limit < MIN_MINOR_LIMIT)

           _minor_limit = MIN_MINOR_LIMIT;

       }

       /// Finds the next entering edge.

       bool findEnteringEdge() {

         Cost min, curr;

         if (_minor_count < _minor_limit && _candidates.size() > 0) {

           // Minor iteration

           ++_minor_count;

           Edge e;

           min = 0;

           for (int i = 0; i < int(_candidates.size()); ++i) {

             e = _candidates[i];

             if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {

               min = curr;

               _ns._in_edge = e;

             }

           }

           if (min < 0) return true;

         }

         // Major iteration

         _candidates.clear();

         EdgeIt e = _next_edge;

         min = 0;

         for ( ; e != INVALID; ++e) {

           if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {

             _candidates.push_back(e);

             if (curr < min) {

               min = curr;

               _ns._in_edge = e;

             }

             if (int(_candidates.size()) == _list_length) break;

           }

         }

         if (int(_candidates.size()) < _list_length) {

           for (e = EdgeIt(_ns._graph); e != _next_edge; ++e) {

             if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {

               _candidates.push_back(e);

               if (curr < min) {

                 min = curr;

                 _ns._in_edge = e;

               }

               if (int(_candidates.size()) == _list_length) break;

             }

           }

         }

         if (_candidates.size() == 0) return false;

         _minor_count = 1;

         _next_edge = ++e;

         return true;

       }

     }; //class CandidateListPivotRule

   private:

     // State constant for edges at their lower bounds.

     static const int STATE_LOWER =  1;

     // State constant for edges in the spanning tree.

     static const int STATE_TREE  =  0;

     // State constant for edges at their upper bounds.

     static const int STATE_UPPER = -1;

     // Constant for the combined pivot rule.

     static const int COMBINED_PIVOT_MAX_DEG = 5;

   private:

     // The directed graph the algorithm runs on

     SGraph _graph;

     // The original graph

     const Graph &_graph_ref;

     // The original lower bound map

     const LowerMap *_lower;

     // The capacity map

     SCapacityMap _capacity;

     // The cost map

     SCostMap _cost;

     // The supply map

     SSupplyMap _supply;

     bool _valid_supply;

     // Edge map of the current flow

     SCapacityMap _flow;

     // Node map of the current potentials

     SPotentialMap _potential;

     // The depth node map of the spanning tree structure

     IntNodeMap _depth;

     // The parent node map of the spanning tree structure

     NodeNodeMap _parent;

     // The pred_edge node map of the spanning tree structure

     EdgeNodeMap _pred_edge;

     // The thread node map of the spanning tree structure

     NodeNodeMap _thread;

     // The forward node map of the spanning tree structure

     BoolNodeMap _forward;

     // The state edge map

     IntEdgeMap _state;

     // The root node of the starting spanning tree

     Node _root;

     // The reduced cost map

     ReducedCostMap _red_cost;

     // Members for handling the original graph

     FlowMap _flow_result;

     PotentialMap _potential_result;

     NodeRefMap _node_ref;

     EdgeRefMap _edge_ref;

     Edge in_edge;

     Edge _in_edge;

     /// \param _graph The directed graph the algorithm runs on.

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

     /// \param _lower The lower bounds of the edges.

     /// \param lower The lower bounds of the edges.

     /// \param _capacity The capacities (upper bounds) of the edges.

     /// \param capacity The capacities (upper bounds) of the edges.

     /// \param _cost The cost (length) values of the edges.

     /// \param cost The cost (length) values of the edges.

     /// \param _supply The supply values of the nodes (signed).

     /// \param supply The supply values of the nodes (signed).

     NetworkSimplex( const Graph &_graph,

     NetworkSimplex( const Graph &graph,

                     const LowerMap &_lower,

                     const LowerMap &lower,

                     const CapacityMap &_capacity,

                     const CapacityMap &capacity,

                     const CostMap &_cost,

                     const CostMap &cost,

                     const SupplyMap &_supply ) :

                     const SupplyMap &supply ) :

       graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),

       _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),

       supply(graph), flow(graph), flow_result(_graph), potential(graph),

       _cost(_graph), _supply(_graph), _flow(_graph),

       potential_result(_graph), depth(graph), parent(graph),

       _potential(_graph), _depth(_graph), _parent(_graph),

       pred_edge(graph), thread(graph), forward(graph), state(graph),

       _pred_edge(_graph), _thread(_graph), _forward(_graph),

       node_ref(graph_ref), edge_ref(graph_ref),

       _state(_graph), _red_cost(_graph, _cost, _potential),

       red_cost(graph, cost, potential)

       _flow_result(graph), _potential_result(graph),

       _node_ref(graph), _edge_ref(graph)

       for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)

       for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)

         sum += _supply[n];

         sum += supply[n];

       if (!(valid_supply = sum == 0)) return;

       if (!(_valid_supply = sum == 0)) return;

       // Copying graph_ref to graph

       // Copying _graph_ref to _graph

       graph.reserveNode(countNodes(graph_ref) + 1);

       _graph.reserveNode(countNodes(_graph_ref) + 1);

       graph.reserveEdge(countEdges(graph_ref) + countNodes(graph_ref));

       _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));

       copyGraph(graph, graph_ref)

       copyGraph(_graph, _graph_ref)

         .edgeMap(cost, _cost)

         .edgeMap(_cost, cost)

         .nodeRef(node_ref)

         .nodeRef(_node_ref)

         .edgeRef(edge_ref)

         .edgeRef(_edge_ref)

       for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {

       for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {

         capacity[edge_ref[e]] = _capacity[e] - _lower[e];

         _capacity[_edge_ref[e]] = capacity[e] - lower[e];

       }

       }

       for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {

       for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {

         Supply s = _supply[n];

         Supply s = supply[n];

         for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)

         for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)

           s += _lower[e];

           s += lower[e];

         for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)

         for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)

           s -= _lower[e];

           s -= lower[e];

         supply[node_ref[n]] = s;

         _supply[_node_ref[n]] = s;

     /// \param _graph The directed graph the algorithm runs on.

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

     /// \param _capacity The capacities (upper bounds) of the edges.

     /// \param capacity The capacities (upper bounds) of the edges.

     /// \param _cost The cost (length) values of the edges.

     /// \param cost The cost (length) values of the edges.

     /// \param _supply The supply values of the nodes (signed).

     /// \param supply The supply values of the nodes (signed).

     NetworkSimplex( const Graph &_graph,

     NetworkSimplex( const Graph &graph,

                     const CapacityMap &_capacity,

                     const CapacityMap &capacity,

                     const CostMap &_cost,

                     const CostMap &cost,

                     const SupplyMap &_supply ) :

                     const SupplyMap &supply ) :

       graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),

       _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),

       supply(graph), flow(graph), flow_result(_graph), potential(graph),

       _cost(_graph), _supply(_graph), _flow(_graph),

       potential_result(_graph), depth(graph), parent(graph),

       _potential(_graph), _depth(_graph), _parent(_graph),

       pred_edge(graph), thread(graph), forward(graph), state(graph),

       _pred_edge(_graph), _thread(_graph), _forward(_graph),

       node_ref(graph_ref), edge_ref(graph_ref),

       _state(_graph), _red_cost(_graph, _cost, _potential),

       red_cost(graph, cost, potential)

       _flow_result(graph), _potential_result(graph),

       _node_ref(graph), _edge_ref(graph)

       for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)

       for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)

         sum += _supply[n];

         sum += supply[n];

       if (!(valid_supply = sum == 0)) return;

       if (!(_valid_supply = sum == 0)) return;

       // Copying graph_ref to graph

       // Copying _graph_ref to graph

       copyGraph(graph, graph_ref)

       copyGraph(_graph, _graph_ref)

         .edgeMap(capacity, _capacity)

         .edgeMap(_capacity, capacity)

         .edgeMap(cost, _cost)

         .edgeMap(_cost, cost)

         .nodeMap(supply, _supply)

         .nodeMap(_supply, supply)

         .nodeRef(node_ref)

         .nodeRef(_node_ref)

         .edgeRef(edge_ref)

         .edgeRef(_edge_ref)

     /// \param _graph The directed graph the algorithm runs on.

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

     /// \param _lower The lower bounds of the edges.

     /// \param lower The lower bounds of the edges.

     /// \param _capacity The capacities (upper bounds) of the edges.

     /// \param capacity The capacities (upper bounds) of the edges.

     /// \param _cost The cost (length) values of the edges.

     /// \param cost The cost (length) values of the edges.

     /// \param _s The source node.

     /// \param s The source node.

     /// \param _t The target node.

     /// \param t The target node.

     /// \param _flow_value The required amount of flow from node \c _s

     /// \param flow_value The required amount of flow from node \c s

     /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t).

     /// to node \c t (i.e. the supply of \c s and the demand of \c t).

     NetworkSimplex( const Graph &_graph,

     NetworkSimplex( const Graph &graph,

                     const LowerMap &_lower,

                     const LowerMap &lower,

                     const CapacityMap &_capacity,

                     const CapacityMap &capacity,

                     const CostMap &_cost,

                     const CostMap &cost,

                     typename Graph::Node _s,

                     typename Graph::Node s,

                     typename Graph::Node _t,

                     typename Graph::Node t,

                     typename SupplyMap::Value _flow_value ) :

                     typename SupplyMap::Value flow_value ) :

       graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),

       _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),

       supply(graph), flow(graph), flow_result(_graph), potential(graph),

       _cost(_graph), _supply(_graph), _flow(_graph),

       potential_result(_graph), depth(graph), parent(graph),

       _potential(_graph), _depth(_graph), _parent(_graph),

       pred_edge(graph), thread(graph), forward(graph), state(graph),

       _pred_edge(_graph), _thread(_graph), _forward(_graph),

       node_ref(graph_ref), edge_ref(graph_ref),

       _state(_graph), _red_cost(_graph, _cost, _potential),

       red_cost(graph, cost, potential)

       _flow_result(graph), _potential_result(graph),

       _node_ref(graph), _edge_ref(graph)

       // Copying graph_ref to graph

       // Copying _graph_ref to graph

       copyGraph(graph, graph_ref)

       copyGraph(_graph, _graph_ref)

         .edgeMap(cost, _cost)

         .edgeMap(_cost, cost)

         .nodeRef(node_ref)

         .nodeRef(_node_ref)

         .edgeRef(edge_ref)

         .edgeRef(_edge_ref)

       for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {

       for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {

         capacity[edge_ref[e]] = _capacity[e] - _lower[e];

         _capacity[_edge_ref[e]] = capacity[e] - lower[e];

       }

       }

       for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {

       for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {

         Supply s = 0;

         Supply sum = 0;

         if (n == _s) s =  _flow_value;

         if (n == s) sum =  flow_value;

         if (n == _t) s = -_flow_value;

         if (n == t) sum = -flow_value;

         for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)

         for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)

           s += _lower[e];

           sum += lower[e];

         for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)

         for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)

           s -= _lower[e];

           sum -= lower[e];

         supply[node_ref[n]] = s;

         _supply[_node_ref[n]] = sum;

       }

       }

       valid_supply = true;

       _valid_supply = true;

     /// \param _graph The directed graph the algorithm runs on.

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

     /// \param _capacity The capacities (upper bounds) of the edges.

     /// \param capacity The capacities (upper bounds) of the edges.

     /// \param _cost The cost (length) values of the edges.

     /// \param cost The cost (length) values of the edges.

     /// \param _s The source node.

     /// \param s The source node.

     /// \param _t The target node.

     /// \param t The target node.

     /// \param _flow_value The required amount of flow from node \c _s

     /// \param flow_value The required amount of flow from node \c s

     /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t).

     /// to node \c t (i.e. the supply of \c s and the demand of \c t).

     NetworkSimplex( const Graph &_graph,

     NetworkSimplex( const Graph &graph,

                     const CapacityMap &_capacity,

                     const CapacityMap &capacity,

                     const CostMap &_cost,

                     const CostMap &cost,

                     typename Graph::Node _s,

                     typename Graph::Node s,

                     typename Graph::Node _t,

                     typename Graph::Node t,

                     typename SupplyMap::Value _flow_value ) :

                     typename SupplyMap::Value flow_value ) :

       graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),

       _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),

       supply(graph, 0), flow(graph), flow_result(_graph), potential(graph),

       _cost(_graph), _supply(_graph, 0), _flow(_graph),

       potential_result(_graph), depth(graph), parent(graph),

       _potential(_graph), _depth(_graph), _parent(_graph),

       pred_edge(graph), thread(graph), forward(graph), state(graph),

       _pred_edge(_graph), _thread(_graph), _forward(_graph),

       node_ref(graph_ref), edge_ref(graph_ref),

       _state(_graph), _red_cost(_graph, _cost, _potential),

       red_cost(graph, cost, potential)

       _flow_result(graph), _potential_result(graph),

       _node_ref(graph), _edge_ref(graph)

       // Copying graph_ref to graph

       // Copying _graph_ref to graph

       copyGraph(graph, graph_ref)

       copyGraph(_graph, _graph_ref)

         .edgeMap(capacity, _capacity)

         .edgeMap(_capacity, capacity)

         .edgeMap(cost, _cost)

         .edgeMap(_cost, cost)

         .nodeRef(node_ref)

         .nodeRef(_node_ref)

         .edgeRef(edge_ref)

         .edgeRef(_edge_ref)

       supply[node_ref[_s]] =  _flow_value;

       _supply[_node_ref[s]] =  flow_value;

       supply[node_ref[_t]] = -_flow_value;

       _supply[_node_ref[t]] = -flow_value;

       valid_supply = true;

       _valid_supply = true;

     bool run() {

     bool run(PivotRuleEnum pivot_rule = COMBINED_PIVOT) {

       return init() && start();

       return init() && start(pivot_rule);

     }

     }

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

     /// \brief Returns a const reference to the edge map storing the

     ///

     /// found flow.

     /// Returns a const reference to the flow map.

     ///

     /// Returns a const reference to the edge map storing the found flow.

       return flow_result;

       return _flow_result;

     }

     }

     /// \brief Returns a const reference to the potential map (the dual

     /// \brief Returns a const reference to the node map storing the

     /// solution).

     /// found potentials (the dual solution).

     ///

     ///

     /// Returns a const reference to the potential map (the dual

     /// Returns a const reference to the node map storing the found

     /// solution).

     /// potentials (the dual solution).

       return potential_result;

       return _potential_result;

       for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)

       for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)

         c += flow_result[e] * cost[edge_ref[e]];

         c += _flow_result[e] * _cost[_edge_ref[e]];

   protected:

   private:

       if (!valid_supply) return false;

       if (!_valid_supply) return false;

       for (EdgeIt e(graph); e != INVALID; ++e) {

       for (EdgeIt e(_graph); e != INVALID; ++e) {

         flow[e] = 0;

         _flow[e] = 0;

         state[e] = LOWER;

         _state[e] = STATE_LOWER;

       root = graph.addNode();

       _root = _graph.addNode();

       parent[root] = INVALID;

       _parent[_root] = INVALID;

       pred_edge[root] = INVALID;

       _pred_edge[_root] = INVALID;

       depth[root] = 0;

       _depth[_root] = 0;

       supply[root] = 0;

       _supply[_root] = 0;

       potential[root] = 0;

       _potential[_root] = 0;

       Supply sum = 0;

       Node last = _root;

       Node last = root;

       for (NodeIt u(graph); u != INVALID; ++u) {

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

         if (u == root) continue;

         if (u == _root) continue;

         thread[last] = u;

         _thread[last] = u;

         parent[u] = root;

         _parent[u] = _root;

         depth[u] = 1;

         _depth[u] = 1;

         sum += supply[u];

         if (_supply[u] >= 0) {

         if (supply[u] >= 0) {

           e = _graph.addEdge(u, _root);

           e = graph.addEdge(u, root);

           _flow[e] = _supply[u];

           flow[e] = supply[u];

           _forward[u] = true;

           forward[u] = true;

           _potential[u] = -max_cost;

           potential[u] = max_cost;

           e = graph.addEdge(root, u);

           e = _graph.addEdge(_root, u);

           flow[e] = -supply[u];

           _flow[e] = -_supply[u];

           forward[u] = false;

           _forward[u] = false;

           potential[u] = -max_cost;

           _potential[u] = max_cost;

         }

         }

         cost[e] = max_cost;

         _cost[e] = max_cost;

         capacity[e] = std::numeric_limits<Capacity>::max();

         _capacity[e] = std::numeric_limits<Capacity>::max();

         state[e] = TREE;

         _state[e] = STATE_TREE;

         pred_edge[u] = e;

         _pred_edge[u] = e;

       }

       }

       thread[last] = root;

       _thread[last] = _root;

 #ifdef EDGE_BLOCK_PIVOT

       // Initializing block_size for the edge block pivot rule

       int edge_num = countEdges(graph);

       block_size = 2 * int(sqrt(countEdges(graph)));

       if (block_size < MIN_BLOCK_SIZE) block_size = MIN_BLOCK_SIZE;

 #endif

 #ifdef CANDIDATE_LIST_PIVOT

       int edge_num = countEdges(graph);

       minor_count = 0;

       list_length = edge_num / LIST_LENGTH_DIV;

       minor_limit = edge_num / MINOR_LIMIT_DIV;

 #endif

 #ifdef SORTED_LIST_PIVOT

       int edge_num = countEdges(graph);

       list_index = 0;

       list_length = edge_num / LIST_LENGTH_DIV;

 #endif

       return sum == 0;

     }

 #ifdef FIRST_ELIGIBLE_PIVOT

     /// \brief Finds entering edge according to the "First Eligible"

     /// pivot rule.

     bool findEnteringEdge(EdgeIt &next_edge) {

       for (EdgeIt e = next_edge; e != INVALID; ++e) {

         if (state[e] * red_cost[e] < 0) {

           in_edge = e;

           next_edge = ++e;

           return true;

         }

       }

       for (EdgeIt e(graph); e != next_edge; ++e) {

         if (state[e] * red_cost[e] < 0) {

           in_edge = e;

           next_edge = ++e;

           return true;

         }

       }

       return false;

     }

 #endif

 #ifdef BEST_ELIGIBLE_PIVOT

     /// \brief Finds entering edge according to the "Best Eligible"

     /// pivot rule.

     bool findEnteringEdge() {

       Cost min = 0;

       for (EdgeIt e(graph); e != INVALID; ++e) {

         if (state[e] * red_cost[e] < min) {

           min = state[e] * red_cost[e];

           in_edge = e;

         }

       }

       return min < 0;

     }

 #endif

 #ifdef EDGE_BLOCK_PIVOT

     /// \brief Finds entering edge according to the "Edge Block"

     /// pivot rule.

     bool findEnteringEdge(EdgeIt &next_edge) {

       // Performing edge block selection

       Cost curr, min = 0;

       EdgeIt min_edge(graph);

       int cnt = 0;

       for (EdgeIt e = next_edge; e != INVALID; ++e) {

         if ((curr = state[e] * red_cost[e]) < min) {

           min = curr;

           min_edge = e;

         }

         if (++cnt == block_size) {

           if (min < 0) break;

           cnt = 0;

         }

       }

       if (!(min < 0)) {

         for (EdgeIt e(graph); e != next_edge; ++e) {

           if ((curr = state[e] * red_cost[e]) < min) {

             min = curr;

             min_edge = e;

           }

           if (++cnt == block_size) {

             if (min < 0) break;

             cnt = 0;

           }

         }

       }

       in_edge = min_edge;

       if ((next_edge = ++min_edge) == INVALID)

         next_edge = EdgeIt(graph);

       return min < 0;

     }

 #endif

 #ifdef CANDIDATE_LIST_PIVOT

     /// \brief Finds entering edge according to the "Candidate List"

     /// pivot rule.

     bool findEnteringEdge() {

       typedef typename std::vector<Edge>::iterator ListIt;

       if (minor_count >= minor_limit || candidates.size() == 0) {

         // Major iteration

         candidates.clear();

         for (EdgeIt e(graph); e != INVALID; ++e) {

           if (state[e] * red_cost[e] < 0) {

             candidates.push_back(e);

             if (candidates.size() == list_length) break;

           }

         }

         if (candidates.size() == 0) return false;

       }

       // Minor iteration

       ++minor_count;

       Cost min = 0;

       Edge e;

       for (int i = 0; i < candidates.size(); ++i) {

         e = candidates[i];

         if (state[e] * red_cost[e] < min) {

           min = state[e] * red_cost[e];

           in_edge = e;

         }

       }

 #endif

     /// Finds the join node.

 #ifdef SORTED_LIST_PIVOT

     /// \brief Functor class to compare edges during sort of the

     /// candidate list.

     class SortFunc

     {

     private:

       const IntEdgeMap &st;

       const ReducedCostMap &rc;

     public:

       SortFunc(const IntEdgeMap &_st, const ReducedCostMap &_rc) :

         st(_st), rc(_rc) {}

       bool operator()(const Edge &e1, const Edge &e2) {

         return st[e1] * rc[e1] < st[e2] * rc[e2];

       }

     };

     /// \brief Finds entering edge according to the "Sorted List"

     /// pivot rule.

     bool findEnteringEdge() {

       static SortFunc sort_func(state, red_cost);

       // Minor iteration

       while (list_index < candidates.size()) {

         in_edge = candidates[list_index++];

         if (state[in_edge] * red_cost[in_edge] < 0) return true;

       }

       // Major iteration

       candidates.clear();

       Cost curr, min = 0;

       for (EdgeIt e(graph); e != INVALID; ++e) {

         if ((curr = state[e] * red_cost[e]) < min / LOWER_DIV) {

           candidates.push_back(e);

           if (curr < min) min = curr;

           if (candidates.size() == list_length) break;

         }

       }

       if (candidates.size() == 0) return false;

       sort(candidates.begin(), candidates.end(), sort_func);

       in_edge = candidates[0];

       list_index = 1;

       return true;

     }

 #endif

     /// \brief Finds the join node.

       // Finding the join node

       Node u = _graph.source(_in_edge);

       Node u = graph.source(in_edge);

       Node v = _graph.target(_in_edge);

       Node v = graph.target(in_edge);

         if (depth[u] == depth[v]) {

         if (_depth[u] == _depth[v]) {

           u = parent[u];

           u = _parent[u];

           v = parent[v];

           v = _parent[v];

         }

         }

         else if (depth[u] > depth[v]) u = parent[u];

         else if (_depth[u] > _depth[v]) u = _parent[u];

         else v = parent[v];

         else v = _parent[v];

       if (state[in_edge] == LOWER) {

       if (_state[_in_edge] == STATE_LOWER) {

         first = graph.source(in_edge);

         first  = _graph.source(_in_edge);

         second  = graph.target(in_edge);

         second = _graph.target(_in_edge);

         first = graph.target(in_edge);

         first  = _graph.target(_in_edge);

         second  = graph.source(in_edge);

         second = _graph.source(_in_edge);

       }

       }

       delta = capacity[in_edge];

       delta = _capacity[_in_edge];

       for (Node u = first; u != join; u = parent[u]) {

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

         e = pred_edge[u];

         e = _pred_edge[u];

         d = forward[u] ? flow[e] : capacity[e] - flow[e];

         d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e];

       for (Node u = second; u != join; u = parent[u]) {

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

         e = pred_edge[u];

         e = _pred_edge[u];

         d = forward[u] ? capacity[e] - flow[e] : flow[e];

         d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e];

     /// \brief Changes \c flow and \c state edge maps.

     /// Changes \c flow and \c state edge maps.

         Capacity val = state[in_edge] * delta;

         Capacity val = _state[_in_edge] * delta;

         flow[in_edge] += val;

         _flow[_in_edge] += val;

         for (Node u = graph.source(in_edge); u != join; u = parent[u]) {

         for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) {

           flow[pred_edge[u]] += forward[u] ? -val : val;

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

         }

         }

         for (Node u = graph.target(in_edge); u != join; u = parent[u]) {

         for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) {

           flow[pred_edge[u]] += forward[u] ? val : -val;

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

         state[in_edge] = TREE;

         _state[_in_edge] = STATE_TREE;

         state[pred_edge[u_out]] =

         _state[_pred_edge[u_out]] =

           (flow[pred_edge[u_out]] == 0) ? LOWER : UPPER;

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

         state[in_edge] = -state[in_edge];

         _state[_in_edge] = -_state[_in_edge];

       }

       }

     }

     }

     /// \brief Updates \c thread and \c parent node maps.

     /// Updates \c thread and \c parent node maps.

       v_out = parent[u_out];

       v_out = _parent[u_out];

         for (u = join; u != u_in && u != v_in; u = thread[u]) ;

         for (u = join; u != u_in && u != v_in; u = _thread[u]) ;

           while (thread[u] != u_out) u = thread[u];

           while (_thread[u] != u_out) u = _thread[u];

       for (u = u_in; depth[thread[u]] > depth[u_in]; u = thread[u]) ;

       for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;

       right = thread[u];

       right = _thread[u];

       if (thread[v_in] == u_out) {

       if (_thread[v_in] == u_out) {

         for (last = u; depth[last] > depth[u_out]; last = thread[last]) ;

         for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;

         if (last == u_out) last = thread[last];

         if (last == u_out) last = _thread[last];

       }

       }

       else last = thread[v_in];

       else last = _thread[v_in];

       thread[v_in] = stem = u_in;

       _thread[v_in] = stem = u_in;

         thread[u] = new_stem = parent[stem];

         _thread[u] = new_stem = _parent[stem];

         for (u = new_stem; thread[u] != stem; u = thread[u]) ;

         for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;

         thread[u] = right;

         _thread[u] = right;

         parent[stem] = par_stem;

         _parent[stem] = par_stem;

         for (u = stem; depth[thread[u]] > depth[stem]; u = thread[u]) ;

         for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;

         right = thread[u];

         right = _thread[u];

       }

       }

       parent[u_out] = par_stem;

       _parent[u_out] = par_stem;

       thread[u] = last;

       _thread[u] = last;

         if (first != v_in) thread[first] = right;

         if (first != v_in) _thread[first] = right;

         for (u = v_out; thread[u] != u_out; u = thread[u]) ;

         for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;