LEMON/LEMON-official Changeset - r652:5232721b3f14

Repositories » LEMON » LEMON-official

Summary
Changelog
Files
Switch To
- loading...
Options
- Lightweight changelog
- Search
Pull Requests

Changeset - r652:5232721b3f14

« 651:8c3112

653:c7d160 »

r652:5232721b3f14

Wed, 25 Mar 2009 15:58:44

[Not Reviewed]

0 comments (0 inline)

kpeter (Peter Kovacs)
kpeter@inf.elte.hu

Rework the interface of NetworkSimplex (#234) The parameters of the problem can be set with separate functions instead of different constructors.

0 3 0

default

3 files changed with 483 insertions and 552 deletions:

lemon/network_simplex.h

416

318

test/min_cost_flow_test.cc

231

tools/dimacs-solver.cc

↑ Collapse diff ↑

lemon/network_simplex.h

Show inline comments

 		/* -*- mode: C++; indent-tabs-mode: nil; -*-
+		 *
 		 * This file is a part of LEMON, a generic C++ optimization library.
+		 *
 		 * Copyright (C) 2003-2009
 		 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
 		 * (Egervary Research Group on Combinatorial Optimization, EGRES).
+		 *
 		 * Permission to use, modify and distribute this software is granted
 		 * provided that this copyright notice appears in all copies. For
 		 * precise terms see the accompanying LICENSE file.
+		 *
 		 * This software is provided "AS IS" with no warranty of any kind,
 		 * express or implied, and with no claim as to its suitability for any
 		 * purpose.
+		 *
 		 */
 		#ifndef LEMON_NETWORK_SIMPLEX_H
 		#define LEMON_NETWORK_SIMPLEX_H
 		/// \ingroup min_cost_flow
 		///
 		/// \file
-		/// \brief Network simplex algorithm for finding a minimum cost flow.
+		/// \brief Network Simplex algorithm for finding a minimum cost flow.
 		#include <vector>
 		#include <limits>
 		#include <algorithm>
 		#include <lemon/core.h>
 		#include <lemon/math.h>
 		namespace lemon {
 		  /// \addtogroup min_cost_flow
 		  /// @{
-		  /// \brief Implementation of the primal network simplex algorithm
+		  /// \brief Implementation of the primal Network Simplex algorithm
 		  /// for finding a \ref min_cost_flow "minimum cost flow".
 		  ///
-		  /// \ref NetworkSimplex implements the primal network simplex algorithm
+		  /// \ref NetworkSimplex implements the primal Network Simplex algorithm
 		  /// for finding a \ref min_cost_flow "minimum cost flow".
 		  ///
 		  /// \tparam Digraph The digraph type the algorithm runs on.
 		  /// \tparam LowerMap The type of the lower bound map.
 		  /// \tparam CapacityMap The type of the capacity (upper bound) map.
 		  /// \tparam CostMap The type of the cost (length) map.
-		  /// \tparam SupplyMap The type of the supply map.
+		  /// \tparam GR The digraph type the algorithm runs on.
 		  /// \tparam V The value type used in the algorithm.
 		  /// By default it is \c int.
 		  ///
 		  /// \warning
 		  /// - Arc capacities and costs should be \e non-negative \e integers.
 		  /// - Supply values should be \e signed \e integers.
 		  /// - The value types of the maps should be convertible to each other.
-		  /// - \c CostMap::Value must be signed type.
+		  /// \warning \c V must be a signed integer type.
 		  ///
 		  /// \note \ref NetworkSimplex provides five different pivot rule
 		  /// implementations that significantly affect the efficiency of the
 		  /// algorithm.
 		  /// By default "Block Search" pivot rule is used, which proved to be
 		  /// by far the most efficient according to our benchmark tests.
 		  /// However another pivot rule can be selected using \ref run()
 		  /// function with the proper parameter.
 		#ifdef DOXYGEN
 		  template < typename Digraph,
 		             typename LowerMap,
 		             typename CapacityMap,
 		             typename CostMap,
 		             typename SupplyMap >
 		#else
 		  template < typename Digraph,
 		             typename LowerMap = typename Digraph::template ArcMap<int>,
 		             typename CapacityMap = typename Digraph::template ArcMap<int>,
 		             typename CostMap = typename Digraph::template ArcMap<int>,
 		             typename SupplyMap = typename Digraph::template NodeMap<int> >
-		#endif
+		  /// \note %NetworkSimplex provides five different pivot rule
 		  /// implementations. For more information see \ref PivotRule.
 		  template <typename GR, typename V = int>
 		  class NetworkSimplex
+		  {
-		    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
+		  public:
 		    typedef typename CapacityMap::Value Capacity;
 		    typedef typename CostMap::Value Cost;
-		    typedef typename SupplyMap::Value Supply;
+		    /// The value type of the algorithm
 		    typedef V Value;
 		    /// The type of the flow map
 		    typedef typename GR::template ArcMap<Value> FlowMap;
 		    /// The type of the potential map
 		    typedef typename GR::template NodeMap<Value> PotentialMap;
 		  public:
 		    /// \brief Enum type for selecting the pivot rule.
 		    ///
 		    /// Enum type for selecting the pivot rule for the \ref run()
 		    /// function.
 		    ///
 		    /// \ref NetworkSimplex provides five different pivot rule
 		    /// implementations that significantly affect the running time
 		    /// of the algorithm.
 		    /// By default \ref BLOCK_SEARCH "Block Search" is used, which
 		    /// proved to be the most efficient and the most robust on various
 		    /// test inputs according to our benchmark tests.
 		    /// However another pivot rule can be selected using the \ref run()
 		    /// function with the proper parameter.
 		    enum PivotRule {
 		      /// The First Eligible pivot rule.
 		      /// The next eligible arc is selected in a wraparound fashion
 		      /// in every iteration.
 		      FIRST_ELIGIBLE,
 		      /// The Best Eligible pivot rule.
 		      /// The best eligible arc is selected in every iteration.
 		      BEST_ELIGIBLE,
 		      /// The Block Search pivot rule.
 		      /// A specified number of arcs are examined in every iteration
 		      /// in a wraparound fashion and the best eligible arc is selected
 		      /// from this block.
 		      BLOCK_SEARCH,
 		      /// The Candidate List pivot rule.
 		      /// In a major iteration a candidate list is built from eligible arcs
 		      /// in a wraparound fashion and in the following minor iterations
 		      /// the best eligible arc is selected from this list.
 		      CANDIDATE_LIST,
 		      /// The Altering Candidate List pivot rule.
 		      /// It is a modified version of the Candidate List method.
 		      /// It keeps only the several best eligible arcs from the former
 		      /// candidate list and extends this list in every iteration.
 		      ALTERING_LIST
 		    };
 		  private:
 		    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
 		    typedef typename GR::template ArcMap<Value> ValueArcMap;
 		    typedef typename GR::template NodeMap<Value> ValueNodeMap;
 		    typedef std::vector<Arc> ArcVector;
 		    typedef std::vector<Node> NodeVector;
 		    typedef std::vector<int> IntVector;
 		    typedef std::vector<bool> BoolVector;
 		    typedef std::vector<Capacity> CapacityVector;
 		    typedef std::vector<Cost> CostVector;
 		    typedef std::vector<Supply> SupplyVector;
 		  public:
 		    /// The type of the flow map
 		    typedef typename Digraph::template ArcMap<Capacity> FlowMap;
 		    /// The type of the potential map
 		    typedef typename Digraph::template NodeMap<Cost> PotentialMap;
 		  public:
 		    /// Enum type for selecting the pivot rule used by \ref run()
 		    enum PivotRuleEnum {
 		      FIRST_ELIGIBLE_PIVOT,
 		      BEST_ELIGIBLE_PIVOT,
 		      BLOCK_SEARCH_PIVOT,
 		      CANDIDATE_LIST_PIVOT,
 		      ALTERING_LIST_PIVOT
 		    };
-		  private:
+		    typedef std::vector<Value> ValueVector;
 		    // State constants for arcs
 		    enum ArcStateEnum {
 		      STATE_UPPER = -1,
 		      STATE_TREE  =  0,
 		      STATE_LOWER =  1
 		    };
 		  private:
 		    // References for the original data
 		    const Digraph &_graph;
 		    const LowerMap *_orig_lower;
 		    const CapacityMap &_orig_cap;
 		    const CostMap &_orig_cost;
 		    const SupplyMap *_orig_supply;
 		    Node _orig_source;
 		    Node _orig_target;
-		    Capacity _orig_flow_value;
+		    // Data related to the underlying digraph
 		    const GR &_graph;
 		    int _node_num;
 		    int _arc_num;
 		    // Parameters of the problem
 		    ValueArcMap *_plower;
 		    ValueArcMap *_pupper;
 		    ValueArcMap *_pcost;
 		    ValueNodeMap *_psupply;
 		    bool _pstsup;
 		    Node _psource, _ptarget;
 		    Value _pstflow;
 		    // Result maps
 		    FlowMap *_flow_map;
 		    PotentialMap *_potential_map;
 		    bool _local_flow;
 		    bool _local_potential;
 		    // The number of nodes and arcs in the original graph
 		    int _node_num;
 		    int _arc_num;
-		    // Data structures for storing the graph
+		    // Data structures for storing the digraph
 		    IntNodeMap _node_id;
 		    ArcVector _arc_ref;
 		    IntVector _source;
 		    IntVector _target;
 		    // Node and arc maps
 		    CapacityVector _cap;
 		    CostVector _cost;
 		    CostVector _supply;
 		    CapacityVector _flow;
 		    CostVector _pi;
 		    // Node and arc data
 		    ValueVector _cap;
 		    ValueVector _cost;
 		    ValueVector _supply;
 		    ValueVector _flow;
 		    ValueVector _pi;
 		    // Data for storing the spanning tree structure
 		    IntVector _parent;
 		    IntVector _pred;
 		    IntVector _thread;
 		    IntVector _rev_thread;
 		    IntVector _succ_num;
 		    IntVector _last_succ;
 		    IntVector _dirty_revs;
 		    BoolVector _forward;
 		    IntVector _state;
 		    int _root;
 		    // Temporary data used in the current pivot iteration
 		    int in_arc, join, u_in, v_in, u_out, v_out;
 		    int first, second, right, last;
 		    int stem, par_stem, new_stem;
-		    Capacity delta;
+		    Value delta;
 		  private:
 		    /// \brief Implementation of the "First Eligible" pivot rule for the
 		    /// \ref NetworkSimplex "network simplex" algorithm.
 		    ///
 		    /// This class implements the "First Eligible" pivot rule
 		    /// for the \ref NetworkSimplex "network simplex" algorithm.
 		    ///
-		    /// For more information see \ref NetworkSimplex::run().
+		    // Implementation of the First Eligible pivot rule
 		    class FirstEligiblePivotRule
+		    {
 		    private:
 		      // References to the NetworkSimplex class
 		      const IntVector  &_source;
 		      const IntVector  &_target;
-		      const CostVector &_cost;
+		      const ValueVector &_cost;
 		      const IntVector  &_state;
-		      const CostVector &_pi;
+		      const ValueVector &_pi;
 		      int &_in_arc;
 		      int _arc_num;
 		      // Pivot rule data
 		      int _next_arc;
 		    public:
-		      /// Constructor
 		      // Constructor
 		      FirstEligiblePivotRule(NetworkSimplex &ns) :
 		        _source(ns._source), _target(ns._target),
 		        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
 		        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
 		      {}
-		      /// Find next entering arc
 		      // Find next entering arc
 		      bool findEnteringArc() {
-		        Cost c;
+		        Value c;
 		        for (int e = _next_arc; e < _arc_num; ++e) {
 		          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		          if (c < 0) {
 		            _in_arc = e;
 		            _next_arc = e + 1;
 		            return true;
+		          }
+		        }
 		        for (int e = 0; e < _next_arc; ++e) {
 		          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		          if (c < 0) {
 		            _in_arc = e;
 		            _next_arc = e + 1;
 		            return true;
+		          }
+		        }
 		        return false;
+		      }
 		    }; //class FirstEligiblePivotRule
 		    /// \brief Implementation of the "Best Eligible" pivot rule for the
 		    /// \ref NetworkSimplex "network simplex" algorithm.
 		    ///
 		    /// This class implements the "Best Eligible" pivot rule
 		    /// for the \ref NetworkSimplex "network simplex" algorithm.
 		    ///
-		    /// For more information see \ref NetworkSimplex::run().
+		    // Implementation of the Best Eligible pivot rule
 		    class BestEligiblePivotRule
+		    {
 		    private:
 		      // References to the NetworkSimplex class
 		      const IntVector  &_source;
 		      const IntVector  &_target;
-		      const CostVector &_cost;
+		      const ValueVector &_cost;
 		      const IntVector  &_state;
-		      const CostVector &_pi;
+		      const ValueVector &_pi;
 		      int &_in_arc;
 		      int _arc_num;
 		    public:
-		      /// Constructor
 		      // Constructor
 		      BestEligiblePivotRule(NetworkSimplex &ns) :
 		        _source(ns._source), _target(ns._target),
 		        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
 		        _in_arc(ns.in_arc), _arc_num(ns._arc_num)
 		      {}
-		      /// Find next entering arc
 		      // Find next entering arc
 		      bool findEnteringArc() {
-		        Cost c, min = 0;
+		        Value c, min = 0;
 		        for (int e = 0; e < _arc_num; ++e) {
 		          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		          if (c < min) {
 		            min = c;
 		            _in_arc = 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.
 		    ///
-		    /// For more information see \ref NetworkSimplex::run().
+		    // Implementation of the Block Search pivot rule
 		    class BlockSearchPivotRule
+		    {
 		    private:
 		      // References to the NetworkSimplex class
 		      const IntVector  &_source;
 		      const IntVector  &_target;
-		      const CostVector &_cost;
+		      const ValueVector &_cost;
 		      const IntVector  &_state;
-		      const CostVector &_pi;
+		      const ValueVector &_pi;
 		      int &_in_arc;
 		      int _arc_num;
 		      // Pivot rule data
 		      int _block_size;
 		      int _next_arc;
 		    public:
-		      /// Constructor
 		      // Constructor
 		      BlockSearchPivotRule(NetworkSimplex &ns) :
 		        _source(ns._source), _target(ns._target),
 		        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
 		        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
+		      {
 		        // The main parameters of the pivot rule
 		        const double BLOCK_SIZE_FACTOR = 2.0;
 		        const int MIN_BLOCK_SIZE = 10;
 		        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
 		                                MIN_BLOCK_SIZE );
+		      }
-		      /// Find next entering arc
 		      // Find next entering arc
 		      bool findEnteringArc() {
-		        Cost c, min = 0;
+		        Value c, min = 0;
 		        int cnt = _block_size;
 		        int e, min_arc = _next_arc;
 		        for (e = _next_arc; e < _arc_num; ++e) {
 		          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		          if (c < min) {
 		            min = c;
 		            min_arc = e;
+		          }
 		          if (--cnt == 0) {
 		            if (min < 0) break;
 		            cnt = _block_size;
+		          }
+		        }
 		        if (min == 0 || cnt > 0) {
 		          for (e = 0; e < _next_arc; ++e) {
 		            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		            if (c < min) {
 		              min = c;
 		              min_arc = e;
+		            }
 		            if (--cnt == 0) {
 		              if (min < 0) break;
 		              cnt = _block_size;
+		            }
+		          }
+		        }
 		        if (min >= 0) return false;
 		        _in_arc = min_arc;
 		        _next_arc = e;
 		        return true;
+		      }
 		    }; //class BlockSearchPivotRule
 		    /// \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.
 		    ///
-		    /// For more information see \ref NetworkSimplex::run().
+		    // Implementation of the Candidate List pivot rule
 		    class CandidateListPivotRule
+		    {
 		    private:
 		      // References to the NetworkSimplex class
 		      const IntVector  &_source;
 		      const IntVector  &_target;
-		      const CostVector &_cost;
+		      const ValueVector &_cost;
 		      const IntVector  &_state;
-		      const CostVector &_pi;
+		      const ValueVector &_pi;
 		      int &_in_arc;
 		      int _arc_num;
 		      // Pivot rule data
 		      IntVector _candidates;
 		      int _list_length, _minor_limit;
 		      int _curr_length, _minor_count;
 		      int _next_arc;
 		    public:
 		      /// Constructor
 		      CandidateListPivotRule(NetworkSimplex &ns) :
 		        _source(ns._source), _target(ns._target),
 		        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
 		        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
+		      {
 		        // The main parameters of the pivot rule
 		        const double LIST_LENGTH_FACTOR = 1.0;
 		        const int MIN_LIST_LENGTH = 10;
 		        const double MINOR_LIMIT_FACTOR = 0.1;
 		        const int MIN_MINOR_LIMIT = 3;
 		        _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_arc_num)),
 		                                 MIN_LIST_LENGTH );
 		        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
 		                                 MIN_MINOR_LIMIT );
 		        _curr_length = _minor_count = 0;
 		        _candidates.resize(_list_length);
+		      }
 		      /// Find next entering arc
 		      bool findEnteringArc() {
-		        Cost min, c;
+		        Value min, c;
 		        int e, min_arc = _next_arc;
 		        if (_curr_length > 0 && _minor_count < _minor_limit) {
 		          // Minor iteration: select the best eligible arc from the
 		          // current candidate list
 		          ++_minor_count;
 		          min = 0;
 		          for (int i = 0; i < _curr_length; ++i) {
 		            e = _candidates[i];
 		            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		            if (c < min) {
 		              min = c;
 		              min_arc = e;
+		            }
 		            if (c >= 0) {
 		              _candidates[i--] = _candidates[--_curr_length];
+		            }
+		          }
 		          if (min < 0) {
 		            _in_arc = min_arc;
 		            return true;
+		          }
+		        }
 		        // Major iteration: build a new candidate list
 		        min = 0;
 		        _curr_length = 0;
 		        for (e = _next_arc; e < _arc_num; ++e) {
 		          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		          if (c < 0) {
 		            _candidates[_curr_length++] = e;
 		            if (c < min) {
 		              min = c;
 		              min_arc = e;
+		            }
 		            if (_curr_length == _list_length) break;
+		          }
+		        }
 		        if (_curr_length < _list_length) {
 		          for (e = 0; e < _next_arc; ++e) {
 		            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		            if (c < 0) {
 		              _candidates[_curr_length++] = e;
 		              if (c < min) {
 		                min = c;
 		                min_arc = e;
+		              }
 		              if (_curr_length == _list_length) break;
+		            }
+		          }
+		        }
 		        if (_curr_length == 0) return false;
 		        _minor_count = 1;
 		        _in_arc = min_arc;
 		        _next_arc = e;
 		        return true;
+		      }
 		    }; //class CandidateListPivotRule
 		    /// \brief Implementation of the "Altering Candidate List" pivot rule
 		    /// for the \ref NetworkSimplex "network simplex" algorithm.
 		    ///
 		    /// This class implements the "Altering Candidate List" pivot rule
 		    /// for the \ref NetworkSimplex "network simplex" algorithm.
 		    ///
-		    /// For more information see \ref NetworkSimplex::run().
+		    // Implementation of the Altering Candidate List pivot rule
 		    class AlteringListPivotRule
+		    {
 		    private:
 		      // References to the NetworkSimplex class
 		      const IntVector  &_source;
 		      const IntVector  &_target;
-		      const CostVector &_cost;
+		      const ValueVector &_cost;
 		      const IntVector  &_state;
-		      const CostVector &_pi;
+		      const ValueVector &_pi;
 		      int &_in_arc;
 		      int _arc_num;
 		      // Pivot rule data
 		      int _block_size, _head_length, _curr_length;
 		      int _next_arc;
 		      IntVector _candidates;
-		      CostVector _cand_cost;
+		      ValueVector _cand_cost;
 		      // Functor class to compare arcs during sort of the candidate list
 		      class SortFunc
+		      {
 		      private:
-		        const CostVector &_map;
+		        const ValueVector &_map;
 		      public:
-		        SortFunc(const CostVector &map) : _map(map) {}
+		        SortFunc(const ValueVector &map) : _map(map) {}
 		        bool operator()(int left, int right) {
 		          return _map[left] > _map[right];
+		        }
 		      };
 		      SortFunc _sort_func;
 		    public:
-		      /// Constructor
 		      // Constructor
 		      AlteringListPivotRule(NetworkSimplex &ns) :
 		        _source(ns._source), _target(ns._target),
 		        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
 		        _in_arc(ns.in_arc), _arc_num(ns._arc_num),
 		        _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
+		      {
 		        // The main parameters of the pivot rule
 		        const double BLOCK_SIZE_FACTOR = 1.5;
 		        const int MIN_BLOCK_SIZE = 10;
 		        const double HEAD_LENGTH_FACTOR = 0.1;
 		        const int MIN_HEAD_LENGTH = 3;
 		        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
 		                                MIN_BLOCK_SIZE );
 		        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
 		                                 MIN_HEAD_LENGTH );
 		        _candidates.resize(_head_length + _block_size);
 		        _curr_length = 0;
+		      }
-		      /// Find next entering arc
 		      // Find next entering arc
 		      bool findEnteringArc() {
 		        // Check the current candidate list
 		        int e;
 		        for (int i = 0; i < _curr_length; ++i) {
 		          e = _candidates[i];
 		          _cand_cost[e] = _state[e] *
 		            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		          if (_cand_cost[e] >= 0) {
 		            _candidates[i--] = _candidates[--_curr_length];
+		          }
+		        }
 		        // Extend the list
 		        int cnt = _block_size;
 		        int last_arc = 0;
 		        int limit = _head_length;
 		        for (int e = _next_arc; e < _arc_num; ++e) {
 		          _cand_cost[e] = _state[e] *
 		            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		          if (_cand_cost[e] < 0) {
 		            _candidates[_curr_length++] = e;
 		            last_arc = e;
+		          }
 		          if (--cnt == 0) {
 		            if (_curr_length > limit) break;
 		            limit = 0;
 		            cnt = _block_size;
+		          }
+		        }
 		        if (_curr_length <= limit) {
 		          for (int e = 0; e < _next_arc; ++e) {
 		            _cand_cost[e] = _state[e] *
 		              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		            if (_cand_cost[e] < 0) {
 		              _candidates[_curr_length++] = e;
 		              last_arc = e;
+		            }
 		            if (--cnt == 0) {
 		              if (_curr_length > limit) break;
 		              limit = 0;
 		              cnt = _block_size;
+		            }
+		          }
+		        }
 		        if (_curr_length == 0) return false;
 		        _next_arc = last_arc + 1;
 		        // Make heap of the candidate list (approximating a partial sort)
 		        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
 		                   _sort_func );
 		        // Pop the first element of the heap
 		        _in_arc = _candidates[0];
 		        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
 		                  _sort_func );
 		        _curr_length = std::min(_head_length, _curr_length - 1);
 		        return true;
+		      }
 		    }; //class AlteringListPivotRule
 		  public:
-		    /// \brief General constructor (with lower bounds).
+		    /// \brief Constructor.
 		    ///
-		    /// General constructor (with lower bounds).
+		    /// Constructor.
 		    ///
 		    /// \param graph The digraph the algorithm runs on.
 		    /// \param lower The lower bounds of the arcs.
 		    /// \param capacity The capacities (upper bounds) of the arcs.
 		    /// \param cost The cost (length) values of the arcs.
 		    /// \param supply The supply values of the nodes (signed).
 		    NetworkSimplex( const Digraph &graph,
 		                    const LowerMap &lower,
 		                    const CapacityMap &capacity,
 		                    const CostMap &cost,
 		                    const SupplyMap &supply ) :
 		      _graph(graph), _orig_lower(&lower), _orig_cap(capacity),
-		      _orig_cost(cost), _orig_supply(&supply),
+		    NetworkSimplex(const GR& graph) :
 		      _graph(graph),
 		      _plower(NULL), _pupper(NULL), _pcost(NULL),
 		      _psupply(NULL), _pstsup(false),
 		      _flow_map(NULL), _potential_map(NULL),
 		      _local_flow(false), _local_potential(false),
 		      _node_id(graph)
 		    {}
 		    /// \brief General constructor (without lower bounds).
 		    ///
 		    /// General constructor (without lower bounds).
 		    ///
 		    /// \param graph The digraph the algorithm runs on.
 		    /// \param capacity The capacities (upper bounds) of the arcs.
 		    /// \param cost The cost (length) values of the arcs.
 		    /// \param supply The supply values of the nodes (signed).
 		    NetworkSimplex( const Digraph &graph,
 		                    const CapacityMap &capacity,
 		                    const CostMap &cost,
 		                    const SupplyMap &supply ) :
 		      _graph(graph), _orig_lower(NULL), _orig_cap(capacity),
 		      _orig_cost(cost), _orig_supply(&supply),
 		      _flow_map(NULL), _potential_map(NULL),
 		      _local_flow(false), _local_potential(false),
 		      _node_id(graph)
 		    {}
 		    /// \brief Simple constructor (with lower bounds).
 		    ///
 		    /// Simple constructor (with lower bounds).
 		    ///
 		    /// \param graph The digraph the algorithm runs on.
 		    /// \param lower The lower bounds of the arcs.
 		    /// \param capacity The capacities (upper bounds) of the arcs.
 		    /// \param cost The cost (length) values of the arcs.
 		    /// \param s The source node.
 		    /// \param t The target node.
 		    /// \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).
 		    NetworkSimplex( const Digraph &graph,
 		                    const LowerMap &lower,
 		                    const CapacityMap &capacity,
 		                    const CostMap &cost,
 		                    Node s, Node t,
 		                    Capacity flow_value ) :
 		      _graph(graph), _orig_lower(&lower), _orig_cap(capacity),
 		      _orig_cost(cost), _orig_supply(NULL),
 		      _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
 		      _flow_map(NULL), _potential_map(NULL),
 		      _local_flow(false), _local_potential(false),
 		      _node_id(graph)
 		    {}
 		    /// \brief Simple constructor (without lower bounds).
 		    ///
 		    /// Simple constructor (without lower bounds).
 		    ///
 		    /// \param graph The digraph the algorithm runs on.
 		    /// \param capacity The capacities (upper bounds) of the arcs.
 		    /// \param cost The cost (length) values of the arcs.
 		    /// \param s The source node.
 		    /// \param t The target node.
 		    /// \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).
 		    NetworkSimplex( const Digraph &graph,
 		                    const CapacityMap &capacity,
 		                    const CostMap &cost,
 		                    Node s, Node t,
 		                    Capacity flow_value ) :
 		      _graph(graph), _orig_lower(NULL), _orig_cap(capacity),
 		      _orig_cost(cost), _orig_supply(NULL),
 		      _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
 		      _flow_map(NULL), _potential_map(NULL),
 		      _local_flow(false), _local_potential(false),
 		      _node_id(graph)
 		    {}
+		    {
 		      LEMON_ASSERT(std::numeric_limits<Value>::is_integer &&
 		                   std::numeric_limits<Value>::is_signed,
 		        "The value type of NetworkSimplex must be a signed integer");
+		    }
 		    /// Destructor.
 		    ~NetworkSimplex() {
 		      if (_local_flow) delete _flow_map;
 		      if (_local_potential) delete _potential_map;
+		    }
 		    /// \brief Set the lower bounds on the arcs.
 		    ///
 		    /// This function sets the lower bounds on the arcs.
 		    /// If neither this function nor \ref boundMaps() is used before
 		    /// calling \ref run(), the lower bounds will be set to zero
 		    /// on all arcs.
 		    ///
 		    /// \param map An arc map storing the lower bounds.
 		    /// Its \c Value type must be convertible to the \c Value type
 		    /// of the algorithm.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template <typename LOWER>
 		    NetworkSimplex& lowerMap(const LOWER& map) {
 		      delete _plower;
 		      _plower = new ValueArcMap(_graph);
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
 		        (*_plower)[a] = map[a];
+		      }
 		      return *this;
+		    }
 		    /// \brief Set the upper bounds (capacities) on the arcs.
 		    ///
 		    /// This function sets the upper bounds (capacities) on the arcs.
 		    /// If none of the functions \ref upperMap(), \ref capacityMap()
 		    /// and \ref boundMaps() is used before calling \ref run(),
 		    /// the upper bounds (capacities) will be set to
 		    /// \c std::numeric_limits<Value>::max() on all arcs.
 		    ///
 		    /// \param map An arc map storing the upper bounds.
 		    /// Its \c Value type must be convertible to the \c Value type
 		    /// of the algorithm.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template<typename UPPER>
 		    NetworkSimplex& upperMap(const UPPER& map) {
 		      delete _pupper;
 		      _pupper = new ValueArcMap(_graph);
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
 		        (*_pupper)[a] = map[a];
+		      }
 		      return *this;
+		    }
 		    /// \brief Set the upper bounds (capacities) on the arcs.
 		    ///
 		    /// This function sets the upper bounds (capacities) on the arcs.
 		    /// It is just an alias for \ref upperMap().
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template<typename CAP>
 		    NetworkSimplex& capacityMap(const CAP& map) {
 		      return upperMap(map);
+		    }
 		    /// \brief Set the lower and upper bounds on the arcs.
 		    ///
 		    /// This function sets the lower and upper bounds on the arcs.
 		    /// If neither this function nor \ref lowerMap() is used before
 		    /// calling \ref run(), the lower bounds will be set to zero
 		    /// on all arcs.
 		    /// If none of the functions \ref upperMap(), \ref capacityMap()
 		    /// and \ref boundMaps() is used before calling \ref run(),
 		    /// the upper bounds (capacities) will be set to
 		    /// \c std::numeric_limits<Value>::max() on all arcs.
 		    ///
 		    /// \param lower An arc map storing the lower bounds.
 		    /// \param upper An arc map storing the upper bounds.
 		    ///
 		    /// The \c Value type of the maps must be convertible to the
 		    /// \c Value type of the algorithm.
 		    ///
 		    /// \note This function is just a shortcut of calling \ref lowerMap()
 		    /// and \ref upperMap() separately.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template <typename LOWER, typename UPPER>
 		    NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
 		      return lowerMap(lower).upperMap(upper);
+		    }
 		    /// \brief Set the costs of the arcs.
 		    ///
 		    /// This function sets the costs of the arcs.
 		    /// If it is not used before calling \ref run(), the costs
 		    /// will be set to \c 1 on all arcs.
 		    ///
 		    /// \param map An arc map storing the costs.
 		    /// Its \c Value type must be convertible to the \c Value type
 		    /// of the algorithm.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template<typename COST>
 		    NetworkSimplex& costMap(const COST& map) {
 		      delete _pcost;
 		      _pcost = new ValueArcMap(_graph);
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
 		        (*_pcost)[a] = map[a];
+		      }
 		      return *this;
+		    }
 		    /// \brief Set the supply values of the nodes.
 		    ///
 		    /// This function sets the supply values of the nodes.
 		    /// If neither this function nor \ref stSupply() is used before
 		    /// calling \ref run(), the supply of each node will be set to zero.
 		    /// (It makes sense only if non-zero lower bounds are given.)
 		    ///
 		    /// \param map A node map storing the supply values.
 		    /// Its \c Value type must be convertible to the \c Value type
 		    /// of the algorithm.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template<typename SUP>
 		    NetworkSimplex& supplyMap(const SUP& map) {
 		      delete _psupply;
 		      _pstsup = false;
 		      _psupply = new ValueNodeMap(_graph);
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        (*_psupply)[n] = map[n];
+		      }
 		      return *this;
+		    }
 		    /// \brief Set single source and target nodes and a supply value.
 		    ///
 		    /// This function sets a single source node and a single target node
 		    /// and the required flow value.
 		    /// If neither this function nor \ref supplyMap() is used before
 		    /// calling \ref run(), the supply of each node will be set to zero.
 		    /// (It makes sense only if non-zero lower bounds are given.)
 		    ///
 		    /// \param s The source node.
 		    /// \param t The target node.
 		    /// \param k 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).
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
 		      delete _psupply;
 		      _psupply = NULL;
 		      _pstsup = true;
 		      _psource = s;
 		      _ptarget = t;
 		      _pstflow = k;
 		      return *this;
+		    }
 		    /// \brief Set the flow map.
 		    ///
 		    /// This function sets the flow map.
 		    /// If it is not used before calling \ref run(), an instance will
 		    /// be allocated automatically. The destructor deallocates this
 		    /// automatically allocated map, of course.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    NetworkSimplex& flowMap(FlowMap &map) {
 		      if (_local_flow) {
 		        delete _flow_map;
 		        _local_flow = false;
+		      }
 		      _flow_map = &map;
 		      return *this;
+		    }
 		    /// \brief Set the potential map.
 		    ///
-		    /// This function sets the potential map.
+		    /// This function sets the potential map, which is used for storing
 		    /// the dual solution.
 		    /// If it is not used before calling \ref run(), an instance will
 		    /// be allocated automatically. The destructor deallocates this
 		    /// automatically allocated map, of course.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    NetworkSimplex& potentialMap(PotentialMap &map) {
 		      if (_local_potential) {
 		        delete _potential_map;
 		        _local_potential = false;
+		      }
 		      _potential_map = &map;
 		      return *this;
+		    }
 		    /// \name Execution control
 		    /// The algorithm can be executed using the
-		    /// \ref lemon::NetworkSimplex::run() "run()" function.
+		    /// \name Execution Control
 		    /// The algorithm can be executed using \ref run().
 		    /// @{
 		    /// \brief Run the algorithm.
 		    ///
 		    /// This function runs the algorithm.
 		    /// The paramters can be specified using \ref lowerMap(),
 		    /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
 		    /// \ref costMap(), \ref supplyMap() and \ref stSupply()
 		    /// functions. For example,
 		    /// \code
 		    ///   NetworkSimplex<ListDigraph> ns(graph);
 		    ///   ns.boundMaps(lower, upper).costMap(cost)
 		    ///     .supplyMap(sup).run();
 		    /// \endcode
 		    ///
 		    /// \param pivot_rule The pivot rule that is used during the
 		    /// algorithm.
 		    ///
 		    /// The available pivot rules:
 		    ///
 		    /// - FIRST_ELIGIBLE_PIVOT The next eligible arc is selected in
 		    /// a wraparound fashion in every iteration
 		    /// (\ref FirstEligiblePivotRule).
 		    ///
 		    /// - BEST_ELIGIBLE_PIVOT The best eligible arc is selected in
 		    /// every iteration (\ref BestEligiblePivotRule).
 		    ///
 		    /// - BLOCK_SEARCH_PIVOT A specified number of arcs are examined in
 		    /// every iteration in a wraparound fashion and the best eligible
 		    /// arc is selected from this block (\ref BlockSearchPivotRule).
 		    ///
 		    /// - CANDIDATE_LIST_PIVOT In a major iteration a candidate list is
 		    /// built from eligible arcs in a wraparound fashion and in the
 		    /// following minor iterations the best eligible arc is selected
 		    /// from this list (\ref CandidateListPivotRule).
 		    ///
 		    /// - ALTERING_LIST_PIVOT It is a modified version of the
 		    /// "Candidate List" pivot rule. It keeps only the several best
 		    /// eligible arcs from the former candidate list and extends this
 		    /// list in every iteration (\ref AlteringListPivotRule).
 		    ///
 		    /// According to our comprehensive benchmark tests the "Block Search"
 		    /// pivot rule proved to be the fastest and the most robust on
-		    /// various test inputs. Thus it is the default option.
+		    /// \param pivot_rule The pivot rule that will be used during the
 		    /// algorithm. For more information see \ref PivotRule.
 		    ///
 		    /// \return \c true if a feasible flow can be found.
-		    bool run(PivotRuleEnum pivot_rule = BLOCK_SEARCH_PIVOT) {
+		    bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
 		      return init() && start(pivot_rule);
+		    }
 		    /// @}
 		    /// \name Query Functions
 		    /// The results of the algorithm can be obtained using these
 		    /// functions.\n
 		    /// \ref lemon::NetworkSimplex::run() "run()" must be called before
 		    /// using them.
 		    /// The \ref run() function must be called before using them.
 		    /// @{
 		    /// \brief Return the total cost of the found flow.
 		    ///
 		    /// This function returns the total cost of the found flow.
 		    /// The complexity of the function is \f$ O(e) \f$.
 		    ///
 		    /// \note The return type of the function can be specified as a
 		    /// template parameter. For example,
 		    /// \code
 		    ///   ns.totalCost<double>();
 		    /// \endcode
 		    /// It is useful if the total cost cannot be stored in the \c Value
 		    /// type of the algorithm, which is the default return type of the
 		    /// function.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    template <typename Num>
 		    Num totalCost() const {
 		      Num c = 0;
 		      if (_pcost) {
 		        for (ArcIt e(_graph); e != INVALID; ++e)
 		          c += (*_flow_map)[e] * (*_pcost)[e];
 		      } else {
 		        for (ArcIt e(_graph); e != INVALID; ++e)
 		          c += (*_flow_map)[e];
+		      }
 		      return c;
+		    }
 		#ifndef DOXYGEN
 		    Value totalCost() const {
 		      return totalCost<Value>();
+		    }
 		#endif
 		    /// \brief Return the flow on the given arc.
 		    ///
 		    /// This function returns the flow on the given arc.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    Value flow(const Arc& a) const {
 		      return (*_flow_map)[a];
+		    }
 		    /// \brief Return a const reference to the flow map.
 		    ///
 		    /// This function returns a const reference to an arc map storing
 		    /// the found flow.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    const FlowMap& flowMap() const {
 		      return *_flow_map;
+		    }
 		    /// \brief Return the potential (dual value) of the given node.
 		    ///
 		    /// This function returns the potential (dual value) of the
 		    /// given node.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    Value potential(const Node& n) const {
 		      return (*_potential_map)[n];
+		    }
 		    /// \brief Return a const reference to the potential map
 		    /// (the dual solution).
 		    ///
 		    /// This function returns a const reference to a node map storing
-		    /// the found potentials (the dual solution).
+		    /// the found potentials, which form the dual solution of the
 		    /// \ref min_cost_flow "minimum cost flow" problem.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    const PotentialMap& potentialMap() const {
 		      return *_potential_map;
+		    }
 		    /// \brief Return the flow on the given arc.
 		    ///
 		    /// This function returns the flow on the given arc.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    Capacity flow(const Arc& arc) const {
 		      return (*_flow_map)[arc];
+		    }
 		    /// \brief Return the potential of the given node.
 		    ///
 		    /// This function returns the potential of the given node.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    Cost potential(const Node& node) const {
 		      return (*_potential_map)[node];
+		    }
 		    /// \brief Return the total cost of the found flow.
 		    ///
 		    /// This function returns the total cost of the found flow.
 		    /// The complexity of the function is \f$ O(e) \f$.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    Cost totalCost() const {
 		      Cost c = 0;
 		      for (ArcIt e(_graph); e != INVALID; ++e)
 		        c += (*_flow_map)[e] * _orig_cost[e];
 		      return c;
+		    }
 		    /// @}
 		  private:
 		    // Initialize internal data structures
 		    bool init() {
 		      // Initialize result maps
 		      if (!_flow_map) {
 		        _flow_map = new FlowMap(_graph);
 		        _local_flow = true;
+		      }
 		      if (!_potential_map) {
 		        _potential_map = new PotentialMap(_graph);
 		        _local_potential = true;
+		      }
 		      // Initialize vectors
 		      _node_num = countNodes(_graph);
 		      _arc_num = countArcs(_graph);
 		      int all_node_num = _node_num + 1;
 		      int all_arc_num = _arc_num + _node_num;
 		      if (_node_num == 0) return false;
 		      _arc_ref.resize(_arc_num);
 		      _source.resize(all_arc_num);
 		      _target.resize(all_arc_num);
 		      _cap.resize(all_arc_num);
 		      _cost.resize(all_arc_num);
 		      _supply.resize(all_node_num);
 		      _flow.resize(all_arc_num, 0);
 		      _pi.resize(all_node_num, 0);
 		      _parent.resize(all_node_num);
 		      _pred.resize(all_node_num);
 		      _forward.resize(all_node_num);
 		      _thread.resize(all_node_num);
 		      _rev_thread.resize(all_node_num);
 		      _succ_num.resize(all_node_num);
 		      _last_succ.resize(all_node_num);
 		      _state.resize(all_arc_num, STATE_LOWER);
 		      // Initialize node related data
 		      bool valid_supply = true;
 		      if (_orig_supply) {
 		        Supply sum = 0;
 		      if (!_pstsup && !_psupply) {
 		        _pstsup = true;
 		        _psource = _ptarget = NodeIt(_graph);
 		        _pstflow = 0;
+		      }
 		      if (_psupply) {
 		        Value sum = 0;
 		        int i = 0;
 		        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
 		          _node_id[n] = i;
-		          _supply[i] = (*_orig_supply)[n];
+		          _supply[i] = (*_psupply)[n];
 		          sum += _supply[i];
+		        }
 		        valid_supply = (sum == 0);
 		      } else {
 		        int i = 0;
 		        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
 		          _node_id[n] = i;
 		          _supply[i] = 0;
+		        }
 		        _supply[_node_id[_orig_source]] =  _orig_flow_value;
 		        _supply[_node_id[_orig_target]] = -_orig_flow_value;
 		        _supply[_node_id[_psource]] =  _pstflow;
 		        _supply[_node_id[_ptarget]]   = -_pstflow;
+		      }
 		      if (!valid_supply) return false;
 		      // Set data for the artificial root node
 		      _root = _node_num;
 		      _parent[_root] = -1;
 		      _pred[_root] = -1;
 		      _thread[_root] = 0;
 		      _rev_thread[0] = _root;
 		      _succ_num[_root] = all_node_num;
 		      _last_succ[_root] = _root - 1;
 		      _supply[_root] = 0;
 		      _pi[_root] = 0;
 		      // Store the arcs in a mixed order
 		      int k = std::max(int(sqrt(_arc_num)), 10);
 		      int i = 0;
 		      for (ArcIt e(_graph); e != INVALID; ++e) {
 		        _arc_ref[i] = e;
 		        if ((i += k) >= _arc_num) i = (i % k) + 1;
+		      }
 		      // Initialize arc maps
 		      if (_pupper && _pcost) {
 		      for (int i = 0; i != _arc_num; ++i) {
 		        Arc e = _arc_ref[i];
 		        _source[i] = _node_id[_graph.source(e)];
 		        _target[i] = _node_id[_graph.target(e)];
 		        _cost[i] = _orig_cost[e];
 		        _cap[i] = _orig_cap[e];
 		          _cap[i] = (*_pupper)[e];
 		          _cost[i] = (*_pcost)[e];
+		        }
 		      } else {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          Arc e = _arc_ref[i];
 		          _source[i] = _node_id[_graph.source(e)];
 		          _target[i] = _node_id[_graph.target(e)];
+		        }
 		        if (_pupper) {
 		          for (int i = 0; i != _arc_num; ++i)
 		            _cap[i] = (*_pupper)[_arc_ref[i]];
 		        } else {
 		          Value val = std::numeric_limits<Value>::max();
 		          for (int i = 0; i != _arc_num; ++i)
 		            _cap[i] = val;
+		        }
 		        if (_pcost) {
 		          for (int i = 0; i != _arc_num; ++i)
 		            _cost[i] = (*_pcost)[_arc_ref[i]];
 		        } else {
 		          for (int i = 0; i != _arc_num; ++i)
 		            _cost[i] = 1;
+		        }
+		      }
 		      // Remove non-zero lower bounds
-		      if (_orig_lower) {
+		      if (_plower) {
 		        for (int i = 0; i != _arc_num; ++i) {
-		          Capacity c = (*_orig_lower)[_arc_ref[i]];
+		          Value c = (*_plower)[_arc_ref[i]];
 		          if (c != 0) {
 		            _cap[i] -= c;
 		            _supply[_source[i]] -= c;
 		            _supply[_target[i]] += c;
+		          }
+		        }
+		      }
 		      // Add artificial arcs and initialize the spanning tree data structure
 		      Cost max_cost = std::numeric_limits<Cost>::max() / 4;
 		      Capacity max_cap = std::numeric_limits<Capacity>::max();
 		      Value max_cap = std::numeric_limits<Value>::max();
 		      Value max_cost = std::numeric_limits<Value>::max() / 4;
 		      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
 		        _thread[u] = u + 1;
 		        _rev_thread[u + 1] = u;
 		        _succ_num[u] = 1;
 		        _last_succ[u] = u;
 		        _parent[u] = _root;
 		        _pred[u] = e;
 		        if (_supply[u] >= 0) {
 		          _flow[e] = _supply[u];
 		          _forward[u] = true;
 		          _pi[u] = -max_cost;
 		        } else {
 		          _flow[e] = -_supply[u];
 		          _forward[u] = false;
 		          _pi[u] = max_cost;
+		        }
 		        _cost[e] = max_cost;
 		        _cap[e] = max_cap;
 		        _state[e] = STATE_TREE;
+		      }
 		      return true;
+		    }
 		    // Find the join node
 		    void findJoinNode() {
 		      int u = _source[in_arc];
 		      int v = _target[in_arc];
 		      while (u != v) {
 		        if (_succ_num[u] < _succ_num[v]) {
 		          u = _parent[u];
 		        } else {
 		          v = _parent[v];
+		        }
+		      }
 		      join = u;
+		    }
 		    // Find the leaving arc of the cycle and returns true if the
 		    // leaving arc is not the same as the entering arc
 		    bool findLeavingArc() {
 		      // Initialize first and second nodes according to the direction
 		      // of the cycle
 		      if (_state[in_arc] == STATE_LOWER) {
 		        first  = _source[in_arc];
 		        second = _target[in_arc];
 		      } else {
 		        first  = _target[in_arc];
 		        second = _source[in_arc];
+		      }
 		      delta = _cap[in_arc];
 		      int result = 0;
-		      Capacity d;
+		      Value d;
 		      int e;
 		      // Search the cycle along the path form the first node to the root
 		      for (int u = first; u != join; u = _parent[u]) {
 		        e = _pred[u];
 		        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
 		        if (d < delta) {
 		          delta = d;
 		          u_out = u;
 		          result = 1;
+		        }
+		      }
 		      // Search the cycle along the path form the second node to the root
 		      for (int u = second; u != join; u = _parent[u]) {
 		        e = _pred[u];
 		        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
 		        if (d <= delta) {
 		          delta = d;
 		          u_out = u;
 		          result = 2;
+		        }
+		      }
 		      if (result == 1) {
 		        u_in = first;
 		        v_in = second;
 		      } else {
 		        u_in = second;
 		        v_in = first;
+		      }
 		      return result != 0;
+		    }
 		    // Change _flow and _state vectors
 		    void changeFlow(bool change) {
 		      // Augment along the cycle
 		      if (delta > 0) {
-		        Capacity val = _state[in_arc] * delta;
+		        Value val = _state[in_arc] * delta;
 		        _flow[in_arc] += val;
 		        for (int u = _source[in_arc]; u != join; u = _parent[u]) {
 		          _flow[_pred[u]] += _forward[u] ? -val : val;
+		        }
 		        for (int u = _target[in_arc]; u != join; u = _parent[u]) {
 		          _flow[_pred[u]] += _forward[u] ? val : -val;
+		        }
+		      }
 		      // Update the state of the entering and leaving arcs
 		      if (change) {
 		        _state[in_arc] = STATE_TREE;
 		        _state[_pred[u_out]] =
 		          (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
 		      } else {
 		        _state[in_arc] = -_state[in_arc];
+		      }
+		    }
 		    // Update the tree structure
 		    void updateTreeStructure() {
 		      int u, w;
 		      int old_rev_thread = _rev_thread[u_out];
 		      int old_succ_num = _succ_num[u_out];
 		      int old_last_succ = _last_succ[u_out];
 		      v_out = _parent[u_out];
 		      u = _last_succ[u_in];  // the last successor of u_in
 		      right = _thread[u];    // the node after it
 		      // Handle the case when old_rev_thread equals to v_in
 		      // (it also means that join and v_out coincide)
 		      if (old_rev_thread == v_in) {
 		        last = _thread[_last_succ[u_out]];
 		      } else {
 		        last = _thread[v_in];
+		      }
 		      // Update _thread and _parent along the stem nodes (i.e. the nodes
 		      // between u_in and u_out, whose parent have to be changed)
 		      _thread[v_in] = stem = u_in;
 		      _dirty_revs.clear();
 		      _dirty_revs.push_back(v_in);
 		      par_stem = v_in;
 		      while (stem != u_out) {
 		        // Insert the next stem node into the thread list
 		        new_stem = _parent[stem];
 		        _thread[u] = new_stem;
 		        _dirty_revs.push_back(u);
 		        // Remove the subtree of stem from the thread list
 		        w = _rev_thread[stem];
 		        _thread[w] = right;
 		        _rev_thread[right] = w;
 		        // Change the parent node and shift stem nodes
 		        _parent[stem] = par_stem;
 		        par_stem = stem;
 		        stem = new_stem;
 		        // Update u and right
 		        u = _last_succ[stem] == _last_succ[par_stem] ?
 		          _rev_thread[par_stem] : _last_succ[stem];
 		        right = _thread[u];
+		      }
 		      _parent[u_out] = par_stem;
 		      _thread[u] = last;
 		      _rev_thread[last] = u;
 		      _last_succ[u_out] = u;
 		      // Remove the subtree of u_out from the thread list except for
 		      // the case when old_rev_thread equals to v_in
 		      // (it also means that join and v_out coincide)
 		      if (old_rev_thread != v_in) {
 		        _thread[old_rev_thread] = right;
 		        _rev_thread[right] = old_rev_thread;
+		      }
 		      // Update _rev_thread using the new _thread values
 		      for (int i = 0; i < int(_dirty_revs.size()); ++i) {
 		        u = _dirty_revs[i];
 		        _rev_thread[_thread[u]] = u;
+		      }
 		      // Update _pred, _forward, _last_succ and _succ_num for the
 		      // stem nodes from u_out to u_in
 		      int tmp_sc = 0, tmp_ls = _last_succ[u_out];
 		      u = u_out;
 		      while (u != u_in) {
 		        w = _parent[u];
 		        _pred[u] = _pred[w];
 		        _forward[u] = !_forward[w];
 		        tmp_sc += _succ_num[u] - _succ_num[w];
 		        _succ_num[u] = tmp_sc;
 		        _last_succ[w] = tmp_ls;
 		        u = w;
+		      }
 		      _pred[u_in] = in_arc;
 		      _forward[u_in] = (u_in == _source[in_arc]);
 		      _succ_num[u_in] = old_succ_num;
 		      // Set limits for updating _last_succ form v_in and v_out
 		      // towards the root
 		      int up_limit_in = -1;
 		      int up_limit_out = -1;
 		      if (_last_succ[join] == v_in) {
 		        up_limit_out = join;
 		      } else {
 		        up_limit_in = join;
+		      }
 		      // Update _last_succ from v_in towards the root
 		      for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
 		           u = _parent[u]) {
 		        _last_succ[u] = _last_succ[u_out];
+		      }
 		      // Update _last_succ from v_out towards the root
 		      if (join != old_rev_thread && v_in != old_rev_thread) {
 		        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
 		             u = _parent[u]) {
 		          _last_succ[u] = old_rev_thread;
+		        }
 		      } else {
 		        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
 		             u = _parent[u]) {
 		          _last_succ[u] = _last_succ[u_out];
+		        }
+		      }
 		      // Update _succ_num from v_in to join
 		      for (u = v_in; u != join; u = _parent[u]) {
 		        _succ_num[u] += old_succ_num;
+		      }
 		      // Update _succ_num from v_out to join
 		      for (u = v_out; u != join; u = _parent[u]) {
 		        _succ_num[u] -= old_succ_num;
+		      }
+		    }
 		    // Update potentials
 		    void updatePotential() {
-		      Cost sigma = _forward[u_in] ?
+		      Value sigma = _forward[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) {
 		        // Update in the upper subtree (which contains the root)
 		        int before = _rev_thread[u_in];
 		        int after = _thread[_last_succ[u_in]];
 		        _thread[before] = after;
 		        _pi[_root] -= sigma;
 		        for (int u = _thread[_root]; u != _root; u = _thread[u]) {
 		          _pi[u] -= sigma;
+		        }
 		        _thread[before] = u_in;
 		      } else {
 		        // Update in the lower subtree (which has been moved)
 		        int end = _thread[_last_succ[u_in]];
 		        for (int u = u_in; u != end; u = _thread[u]) {
 		          _pi[u] += sigma;
+		        }
+		      }
+		    }
 		    // Execute the algorithm
-		    bool start(PivotRuleEnum pivot_rule) {
+		    bool start(PivotRule pivot_rule) {
 		      // Select the pivot rule implementation
 		      switch (pivot_rule) {
-		        case FIRST_ELIGIBLE_PIVOT:
+		        case FIRST_ELIGIBLE:
 		          return start<FirstEligiblePivotRule>();
-		        case BEST_ELIGIBLE_PIVOT:
+		        case BEST_ELIGIBLE:
 		          return start<BestEligiblePivotRule>();
-		        case BLOCK_SEARCH_PIVOT:
+		        case BLOCK_SEARCH:
 		          return start<BlockSearchPivotRule>();
-		        case CANDIDATE_LIST_PIVOT:
+		        case CANDIDATE_LIST:
 		          return start<CandidateListPivotRule>();
-		        case ALTERING_LIST_PIVOT:
+		        case ALTERING_LIST:
 		          return start<AlteringListPivotRule>();
+		      }
 		      return false;
+		    }
-		    template<class PivotRuleImplementation>
+		    template <typename PivotRuleImpl>
 		    bool start() {
-		      PivotRuleImplementation pivot(*this);
+		      PivotRuleImpl pivot(*this);
-		      // Execute the network simplex algorithm
+		      // Execute the Network Simplex algorithm
 		      while (pivot.findEnteringArc()) {
 		        findJoinNode();
 		        bool change = findLeavingArc();
 		        changeFlow(change);
 		        if (change) {
 		          updateTreeStructure();
 		          updatePotential();
+		        }
+		      }
 		      // Check if the flow amount equals zero on all the artificial arcs
 		      for (int e = _arc_num; e != _arc_num + _node_num; ++e) {
 		        if (_flow[e] > 0) return false;
+		      }
 		      // Copy flow values to _flow_map
-		      if (_orig_lower) {
+		      if (_plower) {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          Arc e = _arc_ref[i];
-		          _flow_map->set(e, (*_orig_lower)[e] + _flow[i]);
+		          _flow_map->set(e, (*_plower)[e] + _flow[i]);
+		        }
 		      } else {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          _flow_map->set(_arc_ref[i], _flow[i]);
+		        }
+		      }
 		      // Copy potential values to _potential_map
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        _potential_map->set(n, _pi[_node_id[n]]);
+		      }
 		      return true;
+		    }
 		  }; //class NetworkSimplex
 		  ///@}
 		} //namespace lemon
 		#endif //LEMON_NETWORK_SIMPLEX_H

test/min_cost_flow_test.cc

Show inline comments

 		/* -*- mode: C++; indent-tabs-mode: nil; -*-
+		 *
 		 * This file is a part of LEMON, a generic C++ optimization library.
+		 *
 		 * Copyright (C) 2003-2009
 		 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
 		 * (Egervary Research Group on Combinatorial Optimization, EGRES).
+		 *
 		 * Permission to use, modify and distribute this software is granted
 		 * provided that this copyright notice appears in all copies. For
 		 * precise terms see the accompanying LICENSE file.
+		 *
 		 * This software is provided "AS IS" with no warranty of any kind,
 		 * express or implied, and with no claim as to its suitability for any
 		 * purpose.
+		 *
 		 */
 		#include <iostream>
 		#include <fstream>
 		#include <lemon/list_graph.h>
 		#include <lemon/smart_graph.h>
 		#include <lemon/lgf_reader.h>
 		//#include <lemon/cycle_canceling.h>
 		//#include <lemon/capacity_scaling.h>
 		//#include <lemon/cost_scaling.h>
 		#include <lemon/network_simplex.h>
 		//#include <lemon/min_cost_flow.h>
 		//#include <lemon/min_cost_max_flow.h>
 		#include <lemon/concepts/digraph.h>
 		#include <lemon/concept_check.h>
 		#include "test_tools.h"
 		using namespace lemon;
 		char test_lgf[] =
 		  "@nodes\n"
 		  "label  sup1 sup2 sup3\n"
 		  "    1    20   27    0\n"
 		  "    2    -4    0    0\n"
 		  "    3     0    0    0\n"
 		  "    4     0    0    0\n"
 		  "    5     9    0    0\n"
 		  "    6    -6    0    0\n"
 		  "    7     0    0    0\n"
 		  "    8     0    0    0\n"
 		  "    9     3    0    0\n"
 		  "   10    -2    0    0\n"
 		  "   11     0    0    0\n"
 		  "   12   -20  -27    0\n"
 		  "\n"
 		  "@arcs\n"
 		  "       cost  cap low1 low2\n"
 		  " 1  2    70   11    0    8\n"
 		  " 1  3   150    3    0    1\n"
 		  " 1  4    80   15    0    2\n"
 		  " 2  8    80   12    0    0\n"
 		  " 3  5   140    5    0    3\n"
 		  " 4  6    60   10    0    1\n"
 		  " 4  7    80    2    0    0\n"
 		  " 4  8   110    3    0    0\n"
 		  " 5  7    60   14    0    0\n"
 		  " 5 11   120   12    0    0\n"
 		  " 6  3     0    3    0    0\n"
 		  " 6  9   140    4    0    0\n"
 		  " 6 10    90    8    0    0\n"
 		  " 7  1    30    5    0    0\n"
 		  " 8 12    60   16    0    4\n"
 		  " 9 12    50    6    0    0\n"
 		  "10 12    70   13    0    5\n"
 		  "10  2   100    7    0    0\n"
 		  "10  7    60   10    0    0\n"
 		  "11 10    20   14    0    6\n"
 		  "12 11    30   10    0    0\n"
 		  "\n"
 		  "@attributes\n"
 		  "source 1\n"
 		  "target 12\n";
 		// Check the interface of an MCF algorithm
 		template <typename GR, typename Value>
 		class McfClassConcept
+		{
 		public:
 		  template <typename MCF>
 		  struct Constraints {
 		    void constraints() {
 		      checkConcept<concepts::Digraph, GR>();
 		      MCF mcf_test1(g, lower, upper, cost, sup);
 		      MCF mcf_test2(g, upper, cost, sup);
 		      MCF mcf_test3(g, lower, upper, cost, n, n, k);
 		      MCF mcf_test4(g, upper, cost, n, n, k);
 		      MCF mcf(g);
 		      // TODO: This part should be enabled and the next part
 		      // should be removed if map copying is supported
 		/*
 		      flow = mcf_test1.flowMap();
-		      mcf_test1.flowMap(flow);
+		      b = mcf.lowerMap(lower)
 		             .upperMap(upper)
 		             .capacityMap(upper)
 		             .boundMaps(lower, upper)
 		             .costMap(cost)
 		             .supplyMap(sup)
 		             .stSupply(n, n, k)
 		             .run();
 		      pot = mcf_test1.potentialMap();
 		      mcf_test1.potentialMap(pot);
 		*/
 		/**/
 		      const typename MCF::FlowMap &fm =
 		        mcf_test1.flowMap();
 		      mcf_test1.flowMap(flow);
 		      const typename MCF::PotentialMap &pm =
 		        mcf_test1.potentialMap();
 		      mcf_test1.potentialMap(pot);
 		      const typename MCF::FlowMap &fm = mcf.flowMap();
 		      const typename MCF::PotentialMap &pm = mcf.potentialMap();
 		      v = mcf.totalCost();
 		      double x = mcf.template totalCost<double>();
 		      v = mcf.flow(a);
 		      v = mcf.potential(n);
 		      mcf.flowMap(flow);
 		      mcf.potentialMap(pot);
 		      ignore_unused_variable_warning(fm);
 		      ignore_unused_variable_warning(pm);
 		/**/
 		      mcf_test1.run();
 		      v = mcf_test1.totalCost();
 		      v = mcf_test1.flow(a);
-		      v = mcf_test1.potential(n);
+		      ignore_unused_variable_warning(x);
+		    }
 		    typedef typename GR::Node Node;
 		    typedef typename GR::Arc Arc;
 		    typedef concepts::ReadMap<Node, Value> NM;
 		    typedef concepts::ReadMap<Arc, Value> AM;
 		    const GR &g;
 		    const AM &lower;
 		    const AM &upper;
 		    const AM &cost;
 		    const NM &sup;
 		    const Node &n;
 		    const Arc &a;
 		    const Value &k;
 		    Value v;
 		    bool b;
 		    typename MCF::FlowMap &flow;
 		    typename MCF::PotentialMap &pot;
 		  };
 		};
 		// Check the feasibility of the given flow (primal soluiton)
 		template < typename GR, typename LM, typename UM,
 		           typename SM, typename FM >
 		bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
 		                const SM& supply, const FM& flow )
+		{
 		  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
 		  for (ArcIt e(gr); e != INVALID; ++e) {
 		    if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
+		  }
 		  for (NodeIt n(gr); n != INVALID; ++n) {
 		    typename SM::Value sum = 0;
 		    for (OutArcIt e(gr, n); e != INVALID; ++e)
 		      sum += flow[e];
 		    for (InArcIt e(gr, n); e != INVALID; ++e)
 		      sum -= flow[e];
 		    if (sum != supply[n]) return false;
+		  }
 		  return true;
+		}
 		// Check the feasibility of the given potentials (dual soluiton)
 		// using the Complementary Slackness optimality condition
+		// using the "Complementary Slackness" optimality condition
 		template < typename GR, typename LM, typename UM,
 		           typename CM, typename FM, typename PM >
 		bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
 		                     const CM& cost, const FM& flow, const PM& pi )
+		{
 		  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
 		  bool opt = true;
 		  for (ArcIt e(gr); opt && e != INVALID; ++e) {
 		    typename CM::Value red_cost =
 		      cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
 		    opt = red_cost == 0 ||
 		          (red_cost > 0 && flow[e] == lower[e]) ||
 		          (red_cost < 0 && flow[e] == upper[e]);
+		  }
 		  return opt;
+		}
 		// Run a minimum cost flow algorithm and check the results
 		template < typename MCF, typename GR,
 		           typename LM, typename UM,
 		           typename CM, typename SM >
 		void checkMcf( const MCF& mcf, bool mcf_result,
 		               const GR& gr, const LM& lower, const UM& upper,
 		               const CM& cost, const SM& supply,
 		               bool result, typename CM::Value total,
 		               const std::string &test_id = "" )
+		{
 		  check(mcf_result == result, "Wrong result " + test_id);
 		  if (result) {
 		    check(checkFlow(gr, lower, upper, supply, mcf.flowMap()),
 		          "The flow is not feasible " + test_id);
 		    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
 		    check(checkPotential(gr, lower, upper, cost, mcf.flowMap(),
 		                         mcf.potentialMap()),
 		          "Wrong potentials " + test_id);
+		  }
+		}
 		int main()
+		{
 		  // Check the interfaces
+		  {
 		    typedef int Value;
 		    // This typedef should be enabled if the standard maps are
 		    // reference maps in the graph concepts
 		    // TODO: This typedef should be enabled if the standard maps are
 		    // reference maps in the graph concepts (See #190).
 		/**/
 		    //typedef concepts::Digraph GR;
 		    typedef ListDigraph GR;
 		    typedef concepts::ReadMap<GR::Node, Value> NM;
 		    typedef concepts::ReadMap<GR::Arc, Value> AM;
 		    //checkConcept< McfClassConcept<GR, Value>,
 		    //              CycleCanceling<GR, AM, AM, AM, NM> >();
 		    //checkConcept< McfClassConcept<GR, Value>,
 		    //              CapacityScaling<GR, AM, AM, AM, NM> >();
 		    //checkConcept< McfClassConcept<GR, Value>,
-		    //              CostScaling<GR, AM, AM, AM, NM> >();
+		/**/
 		    checkConcept< McfClassConcept<GR, Value>,
 		                  NetworkSimplex<GR, AM, AM, AM, NM> >();
 		    //checkConcept< MinCostFlow<GR, Value>,
-		    //              NetworkSimplex<GR, AM, AM, AM, NM> >();
+		                  NetworkSimplex<GR, Value> >();
+		  }
 		  // Run various MCF tests
 		  typedef ListDigraph Digraph;
 		  DIGRAPH_TYPEDEFS(ListDigraph);
 		  // Read the test digraph
 		  Digraph gr;
 		  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
 		  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr);
 		  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
 		  Node v, w;
 		  std::istringstream input(test_lgf);
 		  DigraphReader<Digraph>(gr, input)
 		    .arcMap("cost", c)
 		    .arcMap("cap", u)
 		    .arcMap("low1", l1)
 		    .arcMap("low2", l2)
 		    .nodeMap("sup1", s1)
 		    .nodeMap("sup2", s2)
 		    .nodeMap("sup3", s3)
 		    .node("source", v)
 		    .node("target", w)
 		    .run();
 		/*
 		  // A. Test CapacityScaling with scaling
 		  // A. Test NetworkSimplex with the default pivot rule
+		  {
 		    CapacityScaling<Digraph> mcf1(gr, u, c, s1);
 		    CapacityScaling<Digraph> mcf2(gr, u, c, v, w, 27);
 		    CapacityScaling<Digraph> mcf3(gr, u, c, s3);
 		    CapacityScaling<Digraph> mcf4(gr, l2, u, c, s1);
 		    CapacityScaling<Digraph> mcf5(gr, l2, u, c, v, w, 27);
 		    CapacityScaling<Digraph> mcf6(gr, l2, u, c, s3);
 		    NetworkSimplex<Digraph> mcf1(gr), mcf2(gr), mcf3(gr), mcf4(gr),
 		                            mcf5(gr), mcf6(gr), mcf7(gr), mcf8(gr);
 		    checkMcf(mcf1, mcf1.run(), gr, l1, u, c, s1, true,  5240, "#A1");
 		    checkMcf(mcf2, mcf2.run(), gr, l1, u, c, s2, true,  7620, "#A2");
 		    checkMcf(mcf3, mcf3.run(), gr, l1, u, c, s3, true,     0, "#A3");
 		    checkMcf(mcf4, mcf4.run(), gr, l2, u, c, s1, true,  5970, "#A4");
 		    checkMcf(mcf5, mcf5.run(), gr, l2, u, c, s2, true,  8010, "#A5");
 		    checkMcf(mcf6, mcf6.run(), gr, l2, u, c, s3, false,    0, "#A6");
 		    checkMcf(mcf1, mcf1.upperMap(u).costMap(c).supplyMap(s1).run(),
 		             gr, l1, u, c, s1, true,  5240, "#A1");
 		    checkMcf(mcf2, mcf2.upperMap(u).costMap(c).stSupply(v, w, 27).run(),
 		             gr, l1, u, c, s2, true,  7620, "#A2");
 		    checkMcf(mcf3, mcf3.boundMaps(l2, u).costMap(c).supplyMap(s1).run(),
 		             gr, l2, u, c, s1, true,  5970, "#A3");
 		    checkMcf(mcf4, mcf4.boundMaps(l2, u).costMap(c).stSupply(v, w, 27).run(),
 		             gr, l2, u, c, s2, true,  8010, "#A4");
 		    checkMcf(mcf5, mcf5.supplyMap(s1).run(),
 		             gr, l1, cu, cc, s1, true,  74, "#A5");
 		    checkMcf(mcf6, mcf6.stSupply(v, w, 27).lowerMap(l2).run(),
 		             gr, l2, cu, cc, s2, true,  94, "#A6");
 		    checkMcf(mcf7, mcf7.run(),
 		             gr, l1, cu, cc, s3, true,   0, "#A7");
 		    checkMcf(mcf8, mcf8.boundMaps(l2, u).run(),
 		             gr, l2, u, cc, s3, false,   0, "#A8");
+		  }
-		  // B. Test CapacityScaling without scaling
+		  // B. Test NetworkSimplex with each pivot rule
+		  {
 		    CapacityScaling<Digraph> mcf1(gr, u, c, s1);
 		    CapacityScaling<Digraph> mcf2(gr, u, c, v, w, 27);
 		    CapacityScaling<Digraph> mcf3(gr, u, c, s3);
 		    CapacityScaling<Digraph> mcf4(gr, l2, u, c, s1);
 		    CapacityScaling<Digraph> mcf5(gr, l2, u, c, v, w, 27);
 		    CapacityScaling<Digraph> mcf6(gr, l2, u, c, s3);
 		    NetworkSimplex<Digraph> mcf1(gr), mcf2(gr), mcf3(gr), mcf4(gr), mcf5(gr);
 		    NetworkSimplex<Digraph>::PivotRule pr;
 		    checkMcf(mcf1, mcf1.run(false), gr, l1, u, c, s1, true,  5240, "#B1");
 		    checkMcf(mcf2, mcf2.run(false), gr, l1, u, c, s2, true,  7620, "#B2");
 		    checkMcf(mcf3, mcf3.run(false), gr, l1, u, c, s3, true,     0, "#B3");
 		    checkMcf(mcf4, mcf4.run(false), gr, l2, u, c, s1, true,  5970, "#B4");
 		    checkMcf(mcf5, mcf5.run(false), gr, l2, u, c, s2, true,  8010, "#B5");
 		    checkMcf(mcf6, mcf6.run(false), gr, l2, u, c, s3, false,    0, "#B6");
 		    pr = NetworkSimplex<Digraph>::FIRST_ELIGIBLE;
 		    checkMcf(mcf1, mcf1.boundMaps(l2, u).costMap(c).supplyMap(s1).run(pr),
 		             gr, l2, u, c, s1, true,  5970, "#B1");
 		    pr = NetworkSimplex<Digraph>::BEST_ELIGIBLE;
 		    checkMcf(mcf2, mcf2.boundMaps(l2, u).costMap(c).supplyMap(s1).run(pr),
 		             gr, l2, u, c, s1, true,  5970, "#B2");
 		    pr = NetworkSimplex<Digraph>::BLOCK_SEARCH;
 		    checkMcf(mcf3, mcf3.boundMaps(l2, u).costMap(c).supplyMap(s1).run(pr),
 		             gr, l2, u, c, s1, true,  5970, "#B3");
 		    pr = NetworkSimplex<Digraph>::CANDIDATE_LIST;
 		    checkMcf(mcf4, mcf4.boundMaps(l2, u).costMap(c).supplyMap(s1).run(pr),
 		             gr, l2, u, c, s1, true,  5970, "#B4");
 		    pr = NetworkSimplex<Digraph>::ALTERING_LIST;
 		    checkMcf(mcf5, mcf5.boundMaps(l2, u).costMap(c).supplyMap(s1).run(pr),
 		             gr, l2, u, c, s1, true,  5970, "#B5");
+		  }
 		  // C. Test CostScaling using partial augment-relabel method
+		  {
 		    CostScaling<Digraph> mcf1(gr, u, c, s1);
 		    CostScaling<Digraph> mcf2(gr, u, c, v, w, 27);
 		    CostScaling<Digraph> mcf3(gr, u, c, s3);
 		    CostScaling<Digraph> mcf4(gr, l2, u, c, s1);
 		    CostScaling<Digraph> mcf5(gr, l2, u, c, v, w, 27);
 		    CostScaling<Digraph> mcf6(gr, l2, u, c, s3);
 		    checkMcf(mcf1, mcf1.run(), gr, l1, u, c, s1, true,  5240, "#C1");
 		    checkMcf(mcf2, mcf2.run(), gr, l1, u, c, s2, true,  7620, "#C2");
 		    checkMcf(mcf3, mcf3.run(), gr, l1, u, c, s3, true,     0, "#C3");
 		    checkMcf(mcf4, mcf4.run(), gr, l2, u, c, s1, true,  5970, "#C4");
 		    checkMcf(mcf5, mcf5.run(), gr, l2, u, c, s2, true,  8010, "#C5");
 		    checkMcf(mcf6, mcf6.run(), gr, l2, u, c, s3, false,    0, "#C6");
+		  }
 		  // D. Test CostScaling using push-relabel method
+		  {
 		    CostScaling<Digraph> mcf1(gr, u, c, s1);
 		    CostScaling<Digraph> mcf2(gr, u, c, v, w, 27);
 		    CostScaling<Digraph> mcf3(gr, u, c, s3);
 		    CostScaling<Digraph> mcf4(gr, l2, u, c, s1);
 		    CostScaling<Digraph> mcf5(gr, l2, u, c, v, w, 27);
 		    CostScaling<Digraph> mcf6(gr, l2, u, c, s3);
 		    checkMcf(mcf1, mcf1.run(false), gr, l1, u, c, s1, true,  5240, "#D1");
 		    checkMcf(mcf2, mcf2.run(false), gr, l1, u, c, s2, true,  7620, "#D2");
 		    checkMcf(mcf3, mcf3.run(false), gr, l1, u, c, s3, true,     0, "#D3");
 		    checkMcf(mcf4, mcf4.run(false), gr, l2, u, c, s1, true,  5970, "#D4");
 		    checkMcf(mcf5, mcf5.run(false), gr, l2, u, c, s2, true,  8010, "#D5");
 		    checkMcf(mcf6, mcf6.run(false), gr, l2, u, c, s3, false,    0, "#D6");
+		  }
 		*/
 		  // E. Test NetworkSimplex with FIRST_ELIGIBLE_PIVOT
+		  {
 		    NetworkSimplex<Digraph>::PivotRuleEnum pr =
 		      NetworkSimplex<Digraph>::FIRST_ELIGIBLE_PIVOT;
 		    NetworkSimplex<Digraph> mcf1(gr, u, c, s1);
 		    NetworkSimplex<Digraph> mcf2(gr, u, c, v, w, 27);
 		    NetworkSimplex<Digraph> mcf3(gr, u, c, s3);
 		    NetworkSimplex<Digraph> mcf4(gr, l2, u, c, s1);
 		    NetworkSimplex<Digraph> mcf5(gr, l2, u, c, v, w, 27);
 		    NetworkSimplex<Digraph> mcf6(gr, l2, u, c, s3);
 		    checkMcf(mcf1, mcf1.run(pr), gr, l1, u, c, s1, true,  5240, "#E1");
 		    checkMcf(mcf2, mcf2.run(pr), gr, l1, u, c, s2, true,  7620, "#E2");
 		    checkMcf(mcf3, mcf3.run(pr), gr, l1, u, c, s3, true,     0, "#E3");
 		    checkMcf(mcf4, mcf4.run(pr), gr, l2, u, c, s1, true,  5970, "#E4");
 		    checkMcf(mcf5, mcf5.run(pr), gr, l2, u, c, s2, true,  8010, "#E5");
 		    checkMcf(mcf6, mcf6.run(pr), gr, l2, u, c, s3, false,    0, "#E6");
+		  }
 		  // F. Test NetworkSimplex with BEST_ELIGIBLE_PIVOT
+		  {
 		    NetworkSimplex<Digraph>::PivotRuleEnum pr =
 		      NetworkSimplex<Digraph>::BEST_ELIGIBLE_PIVOT;
 		    NetworkSimplex<Digraph> mcf1(gr, u, c, s1);
 		    NetworkSimplex<Digraph> mcf2(gr, u, c, v, w, 27);
 		    NetworkSimplex<Digraph> mcf3(gr, u, c, s3);
 		    NetworkSimplex<Digraph> mcf4(gr, l2, u, c, s1);
 		    NetworkSimplex<Digraph> mcf5(gr, l2, u, c, v, w, 27);
 		    NetworkSimplex<Digraph> mcf6(gr, l2, u, c, s3);
 		    checkMcf(mcf1, mcf1.run(pr), gr, l1, u, c, s1, true,  5240, "#F1");
 		    checkMcf(mcf2, mcf2.run(pr), gr, l1, u, c, s2, true,  7620, "#F2");
 		    checkMcf(mcf3, mcf3.run(pr), gr, l1, u, c, s3, true,     0, "#F3");
 		    checkMcf(mcf4, mcf4.run(pr), gr, l2, u, c, s1, true,  5970, "#F4");
 		    checkMcf(mcf5, mcf5.run(pr), gr, l2, u, c, s2, true,  8010, "#F5");
 		    checkMcf(mcf6, mcf6.run(pr), gr, l2, u, c, s3, false,    0, "#F6");
+		  }
 		  // G. Test NetworkSimplex with BLOCK_SEARCH_PIVOT
+		  {
 		    NetworkSimplex<Digraph>::PivotRuleEnum pr =
 		      NetworkSimplex<Digraph>::BLOCK_SEARCH_PIVOT;
 		    NetworkSimplex<Digraph> mcf1(gr, u, c, s1);
 		    NetworkSimplex<Digraph> mcf2(gr, u, c, v, w, 27);
 		    NetworkSimplex<Digraph> mcf3(gr, u, c, s3);
 		    NetworkSimplex<Digraph> mcf4(gr, l2, u, c, s1);
 		    NetworkSimplex<Digraph> mcf5(gr, l2, u, c, v, w, 27);
 		    NetworkSimplex<Digraph> mcf6(gr, l2, u, c, s3);
 		    checkMcf(mcf1, mcf1.run(pr), gr, l1, u, c, s1, true,  5240, "#G1");
 		    checkMcf(mcf2, mcf2.run(pr), gr, l1, u, c, s2, true,  7620, "#G2");
 		    checkMcf(mcf3, mcf3.run(pr), gr, l1, u, c, s3, true,     0, "#G3");
 		    checkMcf(mcf4, mcf4.run(pr), gr, l2, u, c, s1, true,  5970, "#G4");
 		    checkMcf(mcf5, mcf5.run(pr), gr, l2, u, c, s2, true,  8010, "#G5");
 		    checkMcf(mcf6, mcf6.run(pr), gr, l2, u, c, s3, false,    0, "#G6");
+		  }
 		  // H. Test NetworkSimplex with CANDIDATE_LIST_PIVOT
+		  {
 		    NetworkSimplex<Digraph>::PivotRuleEnum pr =
 		      NetworkSimplex<Digraph>::CANDIDATE_LIST_PIVOT;
 		    NetworkSimplex<Digraph> mcf1(gr, u, c, s1);
 		    NetworkSimplex<Digraph> mcf2(gr, u, c, v, w, 27);
 		    NetworkSimplex<Digraph> mcf3(gr, u, c, s3);
 		    NetworkSimplex<Digraph> mcf4(gr, l2, u, c, s1);
 		    NetworkSimplex<Digraph> mcf5(gr, l2, u, c, v, w, 27);
 		    NetworkSimplex<Digraph> mcf6(gr, l2, u, c, s3);
 		    checkMcf(mcf1, mcf1.run(pr), gr, l1, u, c, s1, true,  5240, "#H1");
 		    checkMcf(mcf2, mcf2.run(pr), gr, l1, u, c, s2, true,  7620, "#H2");
 		    checkMcf(mcf3, mcf3.run(pr), gr, l1, u, c, s3, true,     0, "#H3");
 		    checkMcf(mcf4, mcf4.run(pr), gr, l2, u, c, s1, true,  5970, "#H4");
 		    checkMcf(mcf5, mcf5.run(pr), gr, l2, u, c, s2, true,  8010, "#H5");
 		    checkMcf(mcf6, mcf6.run(pr), gr, l2, u, c, s3, false,    0, "#H6");
+		  }
 		  // I. Test NetworkSimplex with ALTERING_LIST_PIVOT
+		  {
 		    NetworkSimplex<Digraph>::PivotRuleEnum pr =
 		      NetworkSimplex<Digraph>::ALTERING_LIST_PIVOT;
 		    NetworkSimplex<Digraph> mcf1(gr, u, c, s1);
 		    NetworkSimplex<Digraph> mcf2(gr, u, c, v, w, 27);
 		    NetworkSimplex<Digraph> mcf3(gr, u, c, s3);
 		    NetworkSimplex<Digraph> mcf4(gr, l2, u, c, s1);
 		    NetworkSimplex<Digraph> mcf5(gr, l2, u, c, v, w, 27);
 		    NetworkSimplex<Digraph> mcf6(gr, l2, u, c, s3);
 		    checkMcf(mcf1, mcf1.run(pr), gr, l1, u, c, s1, true,  5240, "#I1");
 		    checkMcf(mcf2, mcf2.run(pr), gr, l1, u, c, s2, true,  7620, "#I2");
 		    checkMcf(mcf3, mcf3.run(pr), gr, l1, u, c, s3, true,     0, "#I3");
 		    checkMcf(mcf4, mcf4.run(pr), gr, l2, u, c, s1, true,  5970, "#I4");
 		    checkMcf(mcf5, mcf5.run(pr), gr, l2, u, c, s2, true,  8010, "#I5");
 		    checkMcf(mcf6, mcf6.run(pr), gr, l2, u, c, s3, false,    0, "#I6");
+		  }
 		/*
 		  // J. Test MinCostFlow
+		  {
 		    MinCostFlow<Digraph> mcf1(gr, u, c, s1);
 		    MinCostFlow<Digraph> mcf2(gr, u, c, v, w, 27);
 		    MinCostFlow<Digraph> mcf3(gr, u, c, s3);
 		    MinCostFlow<Digraph> mcf4(gr, l2, u, c, s1);
 		    MinCostFlow<Digraph> mcf5(gr, l2, u, c, v, w, 27);
 		    MinCostFlow<Digraph> mcf6(gr, l2, u, c, s3);
 		    checkMcf(mcf1, mcf1.run(), gr, l1, u, c, s1, true,  5240, "#J1");
 		    checkMcf(mcf2, mcf2.run(), gr, l1, u, c, s2, true,  7620, "#J2");
 		    checkMcf(mcf3, mcf3.run(), gr, l1, u, c, s3, true,     0, "#J3");
 		    checkMcf(mcf4, mcf4.run(), gr, l2, u, c, s1, true,  5970, "#J4");
 		    checkMcf(mcf5, mcf5.run(), gr, l2, u, c, s2, true,  8010, "#J5");
 		    checkMcf(mcf6, mcf6.run(), gr, l2, u, c, s3, false,    0, "#J6");
+		  }
 		*/
 		/*
 		  // K. Test MinCostMaxFlow
+		  {
 		    MinCostMaxFlow<Digraph> mcmf(gr, u, c, v, w);
 		    mcmf.run();
 		    checkMcf(mcmf, true, gr, l1, u, c, s3, true, 7620, "#K1");
+		  }
 		*/
 		  return 0;
+		}

tools/dimacs-solver.cc

Show inline comments

 		/* -*- mode: C++; indent-tabs-mode: nil; -*-
+		 *
 		 * This file is a part of LEMON, a generic C++ optimization library.
+		 *
 		 * Copyright (C) 2003-2009
 		 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
 		 * (Egervary Research Group on Combinatorial Optimization, EGRES).
+		 *
 		 * Permission to use, modify and distribute this software is granted
 		 * provided that this copyright notice appears in all copies. For
 		 * precise terms see the accompanying LICENSE file.
+		 *
 		 * This software is provided "AS IS" with no warranty of any kind,
 		 * express or implied, and with no claim as to its suitability for any
 		 * purpose.
+		 *
 		 */
 		///\ingroup tools
 		///\file
 		///\brief DIMACS problem solver.
 		///
 		/// This program solves various problems given in DIMACS format.
 		///
 		/// See
 		/// \verbatim
 		///  dimacs-solver --help
 		/// \endverbatim
 		/// for more info on usage.
 		///
 		#include <iostream>
 		#include <fstream>
 		#include <cstring>
 		#include <lemon/smart_graph.h>
 		#include <lemon/dimacs.h>
 		#include <lemon/lgf_writer.h>
 		#include <lemon/time_measure.h>
 		#include <lemon/arg_parser.h>
 		#include <lemon/error.h>
 		#include <lemon/dijkstra.h>
 		#include <lemon/preflow.h>
 		#include <lemon/max_matching.h>
 		#include <lemon/network_simplex.h>
 		using namespace lemon;
 		typedef SmartDigraph Digraph;
 		DIGRAPH_TYPEDEFS(Digraph);
 		typedef SmartGraph Graph;
 		template<class Value>
 		void solve_sp(ArgParser &ap, std::istream &is, std::ostream &,
 		              DimacsDescriptor &desc)
+		{
 		  bool report = !ap.given("q");
 		  Digraph g;
 		  Node s;
 		  Digraph::ArcMap<Value> len(g);
 		  Timer t;
 		  t.restart();
 		  readDimacsSp(is, g, len, s, desc);
 		  if(report) std::cerr << "Read the file: " << t << '\n';
 		  t.restart();
 		  Dijkstra<Digraph, Digraph::ArcMap<Value> > dij(g,len);
 		  if(report) std::cerr << "Setup Dijkstra class: " << t << '\n';
 		  t.restart();
 		  dij.run(s);
 		  if(report) std::cerr << "Run Dijkstra: " << t << '\n';
+		}
 		template<class Value>
 		void solve_max(ArgParser &ap, std::istream &is, std::ostream &,
 		              DimacsDescriptor &desc)
+		{
 		  bool report = !ap.given("q");
 		  Digraph g;
 		  Node s,t;
 		  Digraph::ArcMap<Value> cap(g);
 		  Timer ti;
 		  ti.restart();
 		  readDimacsMax(is, g, cap, s, t, desc);
 		  if(report) std::cerr << "Read the file: " << ti << '\n';
 		  ti.restart();
 		  Preflow<Digraph, Digraph::ArcMap<Value> > pre(g,cap,s,t);
 		  if(report) std::cerr << "Setup Preflow class: " << ti << '\n';
 		  ti.restart();
 		  pre.run();
 		  if(report) std::cerr << "Run Preflow: " << ti << '\n';
 		  if(report) std::cerr << "\nMax flow value: " << pre.flowValue() << '\n';
+		}
 		template<class Value>
 		void solve_min(ArgParser &ap, std::istream &is, std::ostream &,
 		               DimacsDescriptor &desc)
+		{
 		  bool report = !ap.given("q");
 		  Digraph g;
 		  Digraph::ArcMap<Value> lower(g), cap(g), cost(g);
 		  Digraph::NodeMap<Value> sup(g);
 		  Timer ti;
 		  ti.restart();
 		  readDimacsMin(is, g, lower, cap, cost, sup, desc);
 		  if (report) std::cerr << "Read the file: " << ti << '\n';
 		  ti.restart();
 		  NetworkSimplex< Digraph, Digraph::ArcMap<Value>, Digraph::ArcMap<Value>,
 		                  Digraph::ArcMap<Value>, Digraph::NodeMap<Value> >
-		    ns(g, lower, cap, cost, sup);
+		  NetworkSimplex<Digraph, Value> ns(g);
 		  ns.lowerMap(lower).capacityMap(cap).costMap(cost).supplyMap(sup);
 		  if (report) std::cerr << "Setup NetworkSimplex class: " << ti << '\n';
 		  ti.restart();
 		  ns.run();
 		  if (report) std::cerr << "Run NetworkSimplex: " << ti << '\n';
 		  if (report) std::cerr << "\nMin flow cost: " << ns.totalCost() << '\n';
+		}
 		void solve_mat(ArgParser &ap, std::istream &is, std::ostream &,
 		              DimacsDescriptor &desc)
+		{
 		  bool report = !ap.given("q");
 		  Graph g;
 		  Timer ti;
 		  ti.restart();
 		  readDimacsMat(is, g, desc);
 		  if(report) std::cerr << "Read the file: " << ti << '\n';
 		  ti.restart();
 		  MaxMatching<Graph> mat(g);
 		  if(report) std::cerr << "Setup MaxMatching class: " << ti << '\n';
 		  ti.restart();
 		  mat.run();
 		  if(report) std::cerr << "Run MaxMatching: " << ti << '\n';
 		  if(report) std::cerr << "\nCardinality of max matching: "
 		                       << mat.matchingSize() << '\n';
+		}
 		template<class Value>
 		void solve(ArgParser &ap, std::istream &is, std::ostream &os,
 		           DimacsDescriptor &desc)
+		{
 		  switch(desc.type)
+		    {
 		    case DimacsDescriptor::MIN:
 		      solve_min<Value>(ap,is,os,desc);
 		      break;
 		    case DimacsDescriptor::MAX:
 		      solve_max<Value>(ap,is,os,desc);
 		      break;
 		    case DimacsDescriptor::SP:
 		      solve_sp<Value>(ap,is,os,desc);
 		      break;
 		    case DimacsDescriptor::MAT:
 		      solve_mat(ap,is,os,desc);
 		      break;
 		    default:
 		      break;
+		    }
+		}
 		int main(int argc, const char *argv[]) {
 		  typedef SmartDigraph Digraph;
 		  typedef Digraph::Arc Arc;
 		  std::string inputName;
 		  std::string outputName;
 		  ArgParser ap(argc, argv);
 		  ap.other("[INFILE [OUTFILE]]",
 		           "If either the INFILE or OUTFILE file is missing the standard\n"
 		           "     input/output will be used instead.")
 		    .boolOption("q", "Do not print any report")
 		    .boolOption("int","Use 'int' for capacities, costs etc. (default)")
 		    .optionGroup("datatype","int")
 		#ifdef HAVE_LONG_LONG
 		    .boolOption("long","Use 'long long' for capacities, costs etc.")
 		    .optionGroup("datatype","long")
 		#endif
 		    .boolOption("double","Use 'double' for capacities, costs etc.")
 		    .optionGroup("datatype","double")
 		    .boolOption("ldouble","Use 'long double' for capacities, costs etc.")
 		    .optionGroup("datatype","ldouble")
 		    .onlyOneGroup("datatype")
 		    .run();
 		  std::ifstream input;
 		  std::ofstream output;
 		  switch(ap.files().size())
+		    {
 		    case 2:
 		      output.open(ap.files()[1].c_str());
 		      if (!output) {
 		        throw IoError("Cannot open the file for writing", ap.files()[1]);
+		      }
 		    case 1:
 		      input.open(ap.files()[0].c_str());
 		      if (!input) {
 		        throw IoError("File cannot be found", ap.files()[0]);
+		      }
 		    case 0:
 		      break;
 		    default:
 		      std::cerr << ap.commandName() << ": too many arguments\n";
 		      return 1;
+		    }
 		  std::istream& is = (ap.files().size()<1 ? std::cin : input);
 		  std::ostream& os = (ap.files().size()<2 ? std::cout : output);
 		  DimacsDescriptor desc = dimacsType(is);
 		  if(!ap.given("q"))
+		    {
 		      std::cout << "Problem type: ";
 		      switch(desc.type)
+		        {
 		        case DimacsDescriptor::MIN:
 		          std::cout << "min";
 		          break;
 		        case DimacsDescriptor::MAX:
 		          std::cout << "max";
 		          break;
 		        case DimacsDescriptor::SP:
 		          std::cout << "sp";
 		        case DimacsDescriptor::MAT:
 		          std::cout << "mat";
 		          break;
 		        default:
 		          exit(1);
 		          break;
+		        }
 		      std::cout << "\nNum of nodes: " << desc.nodeNum;
 		      std::cout << "\nNum of arcs:  " << desc.edgeNum;
 		      std::cout << "\n\n";
+		    }
 		  if(ap.given("double"))

0 comments (0 inline)

RhodeCode

Login to your account