COIN-OR::LEMON - Graph Library

Changeset 892:cbf32bf95954 in lemon-1.2


Ignore:
Timestamp:
05/03/10 10:24:52 (9 years ago)
Author:
Alpar Juttner <alpar@…>
Branch:
1.2
Parents:
890:e26ad33d1fbc (diff), 891:5205145fabf6 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Phase:
public
Message:

Merge bugfix #368 to branch 1.2

Files:
2 edited

Legend:

Unmodified
Added
Removed
  • lemon/network_simplex.h

    r877 r892  
    10781078        ART_COST = std::numeric_limits<Cost>::max() / 2 + 1;
    10791079      } else {
    1080         ART_COST = std::numeric_limits<Cost>::min();
     1080        ART_COST = 0;
    10811081        for (int i = 0; i != _arc_num; ++i) {
    10821082          if (_cost[i] > ART_COST) ART_COST = _cost[i];
     
    15901590      if (_sum_supply == 0) {
    15911591        if (_stype == GEQ) {
    1592           Cost max_pot = std::numeric_limits<Cost>::min();
     1592          Cost max_pot = -std::numeric_limits<Cost>::max();
    15931593          for (int i = 0; i != _node_num; ++i) {
    15941594            if (_pi[i] > max_pot) max_pot = _pi[i];
  • lemon/network_simplex.h

    r891 r892  
    33 * This file is a part of LEMON, a generic C++ optimization library.
    44 *
    5  * Copyright (C) 2003-2009
     5 * Copyright (C) 2003-2010
    66 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    77 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     
    4141  ///
    4242  /// \ref NetworkSimplex implements the primal Network Simplex algorithm
    43   /// for finding a \ref min_cost_flow "minimum cost flow".
    44   /// This algorithm is a specialized version of the linear programming
    45   /// simplex method directly for the minimum cost flow problem.
    46   /// It is one of the most efficient solution methods.
     43  /// for finding a \ref min_cost_flow "minimum cost flow"
     44  /// \ref amo93networkflows, \ref dantzig63linearprog,
     45  /// \ref kellyoneill91netsimplex.
     46  /// This algorithm is a highly efficient specialized version of the
     47  /// linear programming simplex method directly for the minimum cost
     48  /// flow problem.
    4749  ///
    48   /// In general this class is the fastest implementation available
    49   /// in LEMON for the minimum cost flow problem.
    50   /// Moreover it supports both directions of the supply/demand inequality
    51   /// constraints. For more information see \ref SupplyType.
     50  /// In general, %NetworkSimplex is the fastest implementation available
     51  /// in LEMON for this problem.
     52  /// Moreover, it supports both directions of the supply/demand inequality
     53  /// constraints. For more information, see \ref SupplyType.
    5254  ///
    5355  /// Most of the parameters of the problem (except for the digraph)
     
    5759  ///
    5860  /// \tparam GR The digraph type the algorithm runs on.
    59   /// \tparam V The value type used for flow amounts, capacity bounds
    60   /// and supply values in the algorithm. By default it is \c int.
    61   /// \tparam C The value type used for costs and potentials in the
    62   /// algorithm. By default it is the same as \c V.
     61  /// \tparam V The number type used for flow amounts, capacity bounds
     62  /// and supply values in the algorithm. By default, it is \c int.
     63  /// \tparam C The number type used for costs and potentials in the
     64  /// algorithm. By default, it is the same as \c V.
    6365  ///
    64   /// \warning Both value types must be signed and all input data must
     66  /// \warning Both number types must be signed and all input data must
    6567  /// be integer.
    6668  ///
    6769  /// \note %NetworkSimplex provides five different pivot rule
    6870  /// implementations, from which the most efficient one is used
    69   /// by default. For more information see \ref PivotRule.
     71  /// by default. For more information, see \ref PivotRule.
    7072  template <typename GR, typename V = int, typename C = V>
    7173  class NetworkSimplex
     
    9698      UNBOUNDED
    9799    };
    98    
     100
    99101    /// \brief Constants for selecting the type of the supply constraints.
    100102    ///
     
    114116      LEQ
    115117    };
    116    
     118
    117119    /// \brief Constants for selecting the pivot rule.
    118120    ///
     
    123125    /// implementations that significantly affect the running time
    124126    /// of the algorithm.
    125     /// By default \ref BLOCK_SEARCH "Block Search" is used, which
     127    /// By default, \ref BLOCK_SEARCH "Block Search" is used, which
    126128    /// proved to be the most efficient and the most robust on various
    127     /// test inputs according to our benchmark tests.
    128     /// However another pivot rule can be selected using the \ref run()
     129    /// test inputs.
     130    /// However, another pivot rule can be selected using the \ref run()
    129131    /// function with the proper parameter.
    130132    enum PivotRule {
    131133
    132       /// The First Eligible pivot rule.
     134      /// The \e First \e Eligible pivot rule.
    133135      /// The next eligible arc is selected in a wraparound fashion
    134136      /// in every iteration.
    135137      FIRST_ELIGIBLE,
    136138
    137       /// The Best Eligible pivot rule.
     139      /// The \e Best \e Eligible pivot rule.
    138140      /// The best eligible arc is selected in every iteration.
    139141      BEST_ELIGIBLE,
    140142
    141       /// The Block Search pivot rule.
     143      /// The \e Block \e Search pivot rule.
    142144      /// A specified number of arcs are examined in every iteration
    143145      /// in a wraparound fashion and the best eligible arc is selected
     
    145147      BLOCK_SEARCH,
    146148
    147       /// The Candidate List pivot rule.
     149      /// The \e Candidate \e List pivot rule.
    148150      /// In a major iteration a candidate list is built from eligible arcs
    149151      /// in a wraparound fashion and in the following minor iterations
     
    151153      CANDIDATE_LIST,
    152154
    153       /// The Altering Candidate List pivot rule.
     155      /// The \e Altering \e Candidate \e List pivot rule.
    154156      /// It is a modified version of the Candidate List method.
    155157      /// It keeps only the several best eligible arcs from the former
     
    157159      ALTERING_LIST
    158160    };
    159    
     161
    160162  private:
    161163
    162164    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
    163165
    164     typedef std::vector<Arc> ArcVector;
    165     typedef std::vector<Node> NodeVector;
    166166    typedef std::vector<int> IntVector;
    167     typedef std::vector<bool> BoolVector;
    168167    typedef std::vector<Value> ValueVector;
    169168    typedef std::vector<Cost> CostVector;
     169    typedef std::vector<char> BoolVector;
     170    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
    170171
    171172    // State constants for arcs
    172     enum ArcStateEnum {
     173    enum ArcState {
    173174      STATE_UPPER = -1,
    174175      STATE_TREE  =  0,
    175176      STATE_LOWER =  1
    176177    };
     178
     179    typedef std::vector<signed char> StateVector;
     180    // Note: vector<signed char> is used instead of vector<ArcState> for
     181    // efficiency reasons
    177182
    178183  private:
     
    195200    IntVector _source;
    196201    IntVector _target;
     202    bool _arc_mixing;
    197203
    198204    // Node and arc data
     
    214220    IntVector _dirty_revs;
    215221    BoolVector _forward;
    216     IntVector _state;
     222    StateVector _state;
    217223    int _root;
    218224
     
    223229    Value delta;
    224230
     231    const Value MAX;
     232
    225233  public:
    226  
     234
    227235    /// \brief Constant for infinite upper bounds (capacities).
    228236    ///
     
    243251      const IntVector  &_target;
    244252      const CostVector &_cost;
    245       const IntVector &_state;
     253      const StateVector &_state;
    246254      const CostVector &_pi;
    247255      int &_in_arc;
     
    264272      bool findEnteringArc() {
    265273        Cost c;
    266         for (int e = _next_arc; e < _search_arc_num; ++e) {
     274        for (int e = _next_arc; e != _search_arc_num; ++e) {
    267275          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    268276          if (c < 0) {
     
    272280          }
    273281        }
    274         for (int e = 0; e < _next_arc; ++e) {
     282        for (int e = 0; e != _next_arc; ++e) {
    275283          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    276284          if (c < 0) {
     
    295303      const IntVector  &_target;
    296304      const CostVector &_cost;
    297       const IntVector &_state;
     305      const StateVector &_state;
    298306      const CostVector &_pi;
    299307      int &_in_arc;
     
    312320      bool findEnteringArc() {
    313321        Cost c, min = 0;
    314         for (int e = 0; e < _search_arc_num; ++e) {
     322        for (int e = 0; e != _search_arc_num; ++e) {
    315323          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    316324          if (c < min) {
     
    334342      const IntVector  &_target;
    335343      const CostVector &_cost;
    336       const IntVector &_state;
     344      const StateVector &_state;
    337345      const CostVector &_pi;
    338346      int &_in_arc;
     
    353361      {
    354362        // The main parameters of the pivot rule
    355         const double BLOCK_SIZE_FACTOR = 0.5;
     363        const double BLOCK_SIZE_FACTOR = 1.0;
    356364        const int MIN_BLOCK_SIZE = 10;
    357365
     
    365373        Cost c, min = 0;
    366374        int cnt = _block_size;
    367         int e, min_arc = _next_arc;
    368         for (e = _next_arc; e < _search_arc_num; ++e) {
     375        int e;
     376        for (e = _next_arc; e != _search_arc_num; ++e) {
    369377          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    370378          if (c < min) {
    371379            min = c;
    372             min_arc = e;
     380            _in_arc = e;
    373381          }
    374382          if (--cnt == 0) {
    375             if (min < 0) break;
     383            if (min < 0) goto search_end;
    376384            cnt = _block_size;
    377385          }
    378386        }
    379         if (min == 0 || cnt > 0) {
    380           for (e = 0; e < _next_arc; ++e) {
    381             c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    382             if (c < min) {
    383               min = c;
    384               min_arc = e;
    385             }
    386             if (--cnt == 0) {
    387               if (min < 0) break;
    388               cnt = _block_size;
    389             }
     387        for (e = 0; e != _next_arc; ++e) {
     388          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     389          if (c < min) {
     390            min = c;
     391            _in_arc = e;
     392          }
     393          if (--cnt == 0) {
     394            if (min < 0) goto search_end;
     395            cnt = _block_size;
    390396          }
    391397        }
    392398        if (min >= 0) return false;
    393         _in_arc = min_arc;
     399
     400      search_end:
    394401        _next_arc = e;
    395402        return true;
     
    408415      const IntVector  &_target;
    409416      const CostVector &_cost;
    410       const IntVector &_state;
     417      const StateVector &_state;
    411418      const CostVector &_pi;
    412419      int &_in_arc;
     
    429436      {
    430437        // The main parameters of the pivot rule
    431         const double LIST_LENGTH_FACTOR = 1.0;
     438        const double LIST_LENGTH_FACTOR = 0.25;
    432439        const int MIN_LIST_LENGTH = 10;
    433440        const double MINOR_LIMIT_FACTOR = 0.1;
     
    446453      bool findEnteringArc() {
    447454        Cost min, c;
    448         int e, min_arc = _next_arc;
     455        int e;
    449456        if (_curr_length > 0 && _minor_count < _minor_limit) {
    450457          // Minor iteration: select the best eligible arc from the
     
    457464            if (c < min) {
    458465              min = c;
    459               min_arc = e;
     466              _in_arc = e;
    460467            }
    461             if (c >= 0) {
     468            else if (c >= 0) {
    462469              _candidates[i--] = _candidates[--_curr_length];
    463470            }
    464471          }
    465           if (min < 0) {
    466             _in_arc = min_arc;
    467             return true;
    468           }
     472          if (min < 0) return true;
    469473        }
    470474
     
    472476        min = 0;
    473477        _curr_length = 0;
    474         for (e = _next_arc; e < _search_arc_num; ++e) {
     478        for (e = _next_arc; e != _search_arc_num; ++e) {
    475479          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    476480          if (c < 0) {
     
    478482            if (c < min) {
    479483              min = c;
    480               min_arc = e;
     484              _in_arc = e;
    481485            }
    482             if (_curr_length == _list_length) break;
    483           }
    484         }
    485         if (_curr_length < _list_length) {
    486           for (e = 0; e < _next_arc; ++e) {
    487             c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    488             if (c < 0) {
    489               _candidates[_curr_length++] = e;
    490               if (c < min) {
    491                 min = c;
    492                 min_arc = e;
    493               }
    494               if (_curr_length == _list_length) break;
     486            if (_curr_length == _list_length) goto search_end;
     487          }
     488        }
     489        for (e = 0; e != _next_arc; ++e) {
     490          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     491          if (c < 0) {
     492            _candidates[_curr_length++] = e;
     493            if (c < min) {
     494              min = c;
     495              _in_arc = e;
    495496            }
     497            if (_curr_length == _list_length) goto search_end;
    496498          }
    497499        }
    498500        if (_curr_length == 0) return false;
     501
     502      search_end:
    499503        _minor_count = 1;
    500         _in_arc = min_arc;
    501504        _next_arc = e;
    502505        return true;
     
    515518      const IntVector  &_target;
    516519      const CostVector &_cost;
    517       const IntVector &_state;
     520      const StateVector &_state;
    518521      const CostVector &_pi;
    519522      int &_in_arc;
     
    550553      {
    551554        // The main parameters of the pivot rule
    552         const double BLOCK_SIZE_FACTOR = 1.5;
     555        const double BLOCK_SIZE_FACTOR = 1.0;
    553556        const int MIN_BLOCK_SIZE = 10;
    554557        const double HEAD_LENGTH_FACTOR = 0.1;
     
    568571        // Check the current candidate list
    569572        int e;
    570         for (int i = 0; i < _curr_length; ++i) {
     573        for (int i = 0; i != _curr_length; ++i) {
    571574          e = _candidates[i];
    572575          _cand_cost[e] = _state[e] *
     
    579582        // Extend the list
    580583        int cnt = _block_size;
    581         int last_arc = 0;
    582584        int limit = _head_length;
    583585
    584         for (int e = _next_arc; e < _search_arc_num; ++e) {
     586        for (e = _next_arc; e != _search_arc_num; ++e) {
    585587          _cand_cost[e] = _state[e] *
    586588            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    587589          if (_cand_cost[e] < 0) {
    588590            _candidates[_curr_length++] = e;
    589             last_arc = e;
    590591          }
    591592          if (--cnt == 0) {
    592             if (_curr_length > limit) break;
     593            if (_curr_length > limit) goto search_end;
    593594            limit = 0;
    594595            cnt = _block_size;
    595596          }
    596597        }
    597         if (_curr_length <= limit) {
    598           for (int e = 0; e < _next_arc; ++e) {
    599             _cand_cost[e] = _state[e] *
    600               (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    601             if (_cand_cost[e] < 0) {
    602               _candidates[_curr_length++] = e;
    603               last_arc = e;
    604             }
    605             if (--cnt == 0) {
    606               if (_curr_length > limit) break;
    607               limit = 0;
    608               cnt = _block_size;
    609             }
     598        for (e = 0; e != _next_arc; ++e) {
     599          _cand_cost[e] = _state[e] *
     600            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     601          if (_cand_cost[e] < 0) {
     602            _candidates[_curr_length++] = e;
     603          }
     604          if (--cnt == 0) {
     605            if (_curr_length > limit) goto search_end;
     606            limit = 0;
     607            cnt = _block_size;
    610608          }
    611609        }
    612610        if (_curr_length == 0) return false;
    613         _next_arc = last_arc + 1;
     611
     612      search_end:
    614613
    615614        // Make heap of the candidate list (approximating a partial sort)
     
    619618        // Pop the first element of the heap
    620619        _in_arc = _candidates[0];
     620        _next_arc = e;
    621621        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
    622622                  _sort_func );
     
    634634    ///
    635635    /// \param graph The digraph the algorithm runs on.
    636     NetworkSimplex(const GR& graph) :
     636    /// \param arc_mixing Indicate if the arcs have to be stored in a
     637    /// mixed order in the internal data structure.
     638    /// In special cases, it could lead to better overall performance,
     639    /// but it is usually slower. Therefore it is disabled by default.
     640    NetworkSimplex(const GR& graph, bool arc_mixing = false) :
    637641      _graph(graph), _node_id(graph), _arc_id(graph),
     642      _arc_mixing(arc_mixing),
     643      MAX(std::numeric_limits<Value>::max()),
    638644      INF(std::numeric_limits<Value>::has_infinity ?
    639           std::numeric_limits<Value>::infinity() :
    640           std::numeric_limits<Value>::max())
     645          std::numeric_limits<Value>::infinity() : MAX)
    641646    {
    642       // Check the value types
     647      // Check the number types
    643648      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
    644649        "The flow type of NetworkSimplex must be signed");
    645650      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
    646651        "The cost type of NetworkSimplex must be signed");
    647        
     652
     653      // Reset data structures
     654      reset();
     655    }
     656
     657    /// \name Parameters
     658    /// The parameters of the algorithm can be specified using these
     659    /// functions.
     660
     661    /// @{
     662
     663    /// \brief Set the lower bounds on the arcs.
     664    ///
     665    /// This function sets the lower bounds on the arcs.
     666    /// If it is not used before calling \ref run(), the lower bounds
     667    /// will be set to zero on all arcs.
     668    ///
     669    /// \param map An arc map storing the lower bounds.
     670    /// Its \c Value type must be convertible to the \c Value type
     671    /// of the algorithm.
     672    ///
     673    /// \return <tt>(*this)</tt>
     674    template <typename LowerMap>
     675    NetworkSimplex& lowerMap(const LowerMap& map) {
     676      _have_lower = true;
     677      for (ArcIt a(_graph); a != INVALID; ++a) {
     678        _lower[_arc_id[a]] = map[a];
     679      }
     680      return *this;
     681    }
     682
     683    /// \brief Set the upper bounds (capacities) on the arcs.
     684    ///
     685    /// This function sets the upper bounds (capacities) on the arcs.
     686    /// If it is not used before calling \ref run(), the upper bounds
     687    /// will be set to \ref INF on all arcs (i.e. the flow value will be
     688    /// unbounded from above).
     689    ///
     690    /// \param map An arc map storing the upper bounds.
     691    /// Its \c Value type must be convertible to the \c Value type
     692    /// of the algorithm.
     693    ///
     694    /// \return <tt>(*this)</tt>
     695    template<typename UpperMap>
     696    NetworkSimplex& upperMap(const UpperMap& map) {
     697      for (ArcIt a(_graph); a != INVALID; ++a) {
     698        _upper[_arc_id[a]] = map[a];
     699      }
     700      return *this;
     701    }
     702
     703    /// \brief Set the costs of the arcs.
     704    ///
     705    /// This function sets the costs of the arcs.
     706    /// If it is not used before calling \ref run(), the costs
     707    /// will be set to \c 1 on all arcs.
     708    ///
     709    /// \param map An arc map storing the costs.
     710    /// Its \c Value type must be convertible to the \c Cost type
     711    /// of the algorithm.
     712    ///
     713    /// \return <tt>(*this)</tt>
     714    template<typename CostMap>
     715    NetworkSimplex& costMap(const CostMap& map) {
     716      for (ArcIt a(_graph); a != INVALID; ++a) {
     717        _cost[_arc_id[a]] = map[a];
     718      }
     719      return *this;
     720    }
     721
     722    /// \brief Set the supply values of the nodes.
     723    ///
     724    /// This function sets the supply values of the nodes.
     725    /// If neither this function nor \ref stSupply() is used before
     726    /// calling \ref run(), the supply of each node will be set to zero.
     727    ///
     728    /// \param map A node map storing the supply values.
     729    /// Its \c Value type must be convertible to the \c Value type
     730    /// of the algorithm.
     731    ///
     732    /// \return <tt>(*this)</tt>
     733    template<typename SupplyMap>
     734    NetworkSimplex& supplyMap(const SupplyMap& map) {
     735      for (NodeIt n(_graph); n != INVALID; ++n) {
     736        _supply[_node_id[n]] = map[n];
     737      }
     738      return *this;
     739    }
     740
     741    /// \brief Set single source and target nodes and a supply value.
     742    ///
     743    /// This function sets a single source node and a single target node
     744    /// and the required flow value.
     745    /// If neither this function nor \ref supplyMap() is used before
     746    /// calling \ref run(), the supply of each node will be set to zero.
     747    ///
     748    /// Using this function has the same effect as using \ref supplyMap()
     749    /// with such a map in which \c k is assigned to \c s, \c -k is
     750    /// assigned to \c t and all other nodes have zero supply value.
     751    ///
     752    /// \param s The source node.
     753    /// \param t The target node.
     754    /// \param k The required amount of flow from node \c s to node \c t
     755    /// (i.e. the supply of \c s and the demand of \c t).
     756    ///
     757    /// \return <tt>(*this)</tt>
     758    NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
     759      for (int i = 0; i != _node_num; ++i) {
     760        _supply[i] = 0;
     761      }
     762      _supply[_node_id[s]] =  k;
     763      _supply[_node_id[t]] = -k;
     764      return *this;
     765    }
     766
     767    /// \brief Set the type of the supply constraints.
     768    ///
     769    /// This function sets the type of the supply/demand constraints.
     770    /// If it is not used before calling \ref run(), the \ref GEQ supply
     771    /// type will be used.
     772    ///
     773    /// For more information, see \ref SupplyType.
     774    ///
     775    /// \return <tt>(*this)</tt>
     776    NetworkSimplex& supplyType(SupplyType supply_type) {
     777      _stype = supply_type;
     778      return *this;
     779    }
     780
     781    /// @}
     782
     783    /// \name Execution Control
     784    /// The algorithm can be executed using \ref run().
     785
     786    /// @{
     787
     788    /// \brief Run the algorithm.
     789    ///
     790    /// This function runs the algorithm.
     791    /// The paramters can be specified using functions \ref lowerMap(),
     792    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
     793    /// \ref supplyType().
     794    /// For example,
     795    /// \code
     796    ///   NetworkSimplex<ListDigraph> ns(graph);
     797    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
     798    ///     .supplyMap(sup).run();
     799    /// \endcode
     800    ///
     801    /// This function can be called more than once. All the given parameters
     802    /// are kept for the next call, unless \ref resetParams() or \ref reset()
     803    /// is used, thus only the modified parameters have to be set again.
     804    /// If the underlying digraph was also modified after the construction
     805    /// of the class (or the last \ref reset() call), then the \ref reset()
     806    /// function must be called.
     807    ///
     808    /// \param pivot_rule The pivot rule that will be used during the
     809    /// algorithm. For more information, see \ref PivotRule.
     810    ///
     811    /// \return \c INFEASIBLE if no feasible flow exists,
     812    /// \n \c OPTIMAL if the problem has optimal solution
     813    /// (i.e. it is feasible and bounded), and the algorithm has found
     814    /// optimal flow and node potentials (primal and dual solutions),
     815    /// \n \c UNBOUNDED if the objective function of the problem is
     816    /// unbounded, i.e. there is a directed cycle having negative total
     817    /// cost and infinite upper bound.
     818    ///
     819    /// \see ProblemType, PivotRule
     820    /// \see resetParams(), reset()
     821    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
     822      if (!init()) return INFEASIBLE;
     823      return start(pivot_rule);
     824    }
     825
     826    /// \brief Reset all the parameters that have been given before.
     827    ///
     828    /// This function resets all the paramaters that have been given
     829    /// before using functions \ref lowerMap(), \ref upperMap(),
     830    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
     831    ///
     832    /// It is useful for multiple \ref run() calls. Basically, all the given
     833    /// parameters are kept for the next \ref run() call, unless
     834    /// \ref resetParams() or \ref reset() is used.
     835    /// If the underlying digraph was also modified after the construction
     836    /// of the class or the last \ref reset() call, then the \ref reset()
     837    /// function must be used, otherwise \ref resetParams() is sufficient.
     838    ///
     839    /// For example,
     840    /// \code
     841    ///   NetworkSimplex<ListDigraph> ns(graph);
     842    ///
     843    ///   // First run
     844    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
     845    ///     .supplyMap(sup).run();
     846    ///
     847    ///   // Run again with modified cost map (resetParams() is not called,
     848    ///   // so only the cost map have to be set again)
     849    ///   cost[e] += 100;
     850    ///   ns.costMap(cost).run();
     851    ///
     852    ///   // Run again from scratch using resetParams()
     853    ///   // (the lower bounds will be set to zero on all arcs)
     854    ///   ns.resetParams();
     855    ///   ns.upperMap(capacity).costMap(cost)
     856    ///     .supplyMap(sup).run();
     857    /// \endcode
     858    ///
     859    /// \return <tt>(*this)</tt>
     860    ///
     861    /// \see reset(), run()
     862    NetworkSimplex& resetParams() {
     863      for (int i = 0; i != _node_num; ++i) {
     864        _supply[i] = 0;
     865      }
     866      for (int i = 0; i != _arc_num; ++i) {
     867        _lower[i] = 0;
     868        _upper[i] = INF;
     869        _cost[i] = 1;
     870      }
     871      _have_lower = false;
     872      _stype = GEQ;
     873      return *this;
     874    }
     875
     876    /// \brief Reset the internal data structures and all the parameters
     877    /// that have been given before.
     878    ///
     879    /// This function resets the internal data structures and all the
     880    /// paramaters that have been given before using functions \ref lowerMap(),
     881    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
     882    /// \ref supplyType().
     883    ///
     884    /// It is useful for multiple \ref run() calls. Basically, all the given
     885    /// parameters are kept for the next \ref run() call, unless
     886    /// \ref resetParams() or \ref reset() is used.
     887    /// If the underlying digraph was also modified after the construction
     888    /// of the class or the last \ref reset() call, then the \ref reset()
     889    /// function must be used, otherwise \ref resetParams() is sufficient.
     890    ///
     891    /// See \ref resetParams() for examples.
     892    ///
     893    /// \return <tt>(*this)</tt>
     894    ///
     895    /// \see resetParams(), run()
     896    NetworkSimplex& reset() {
    648897      // Resize vectors
    649898      _node_num = countNodes(_graph);
     
    672921      _state.resize(max_arc_num);
    673922
    674       // Copy the graph (store the arcs in a mixed order)
     923      // Copy the graph
    675924      int i = 0;
    676925      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
    677926        _node_id[n] = i;
    678927      }
    679       int k = std::max(int(std::sqrt(double(_arc_num))), 10);
    680       i = 0;
    681       for (ArcIt a(_graph); a != INVALID; ++a) {
    682         _arc_id[a] = i;
    683         _source[i] = _node_id[_graph.source(a)];
    684         _target[i] = _node_id[_graph.target(a)];
    685         if ((i += k) >= _arc_num) i = (i % k) + 1;
    686       }
    687      
    688       // Initialize maps
    689       for (int i = 0; i != _node_num; ++i) {
    690         _supply[i] = 0;
    691       }
    692       for (int i = 0; i != _arc_num; ++i) {
    693         _lower[i] = 0;
    694         _upper[i] = INF;
    695         _cost[i] = 1;
    696       }
    697       _have_lower = false;
    698       _stype = GEQ;
    699     }
    700 
    701     /// \name Parameters
    702     /// The parameters of the algorithm can be specified using these
    703     /// functions.
    704 
    705     /// @{
    706 
    707     /// \brief Set the lower bounds on the arcs.
    708     ///
    709     /// This function sets the lower bounds on the arcs.
    710     /// If it is not used before calling \ref run(), the lower bounds
    711     /// will be set to zero on all arcs.
    712     ///
    713     /// \param map An arc map storing the lower bounds.
    714     /// Its \c Value type must be convertible to the \c Value type
    715     /// of the algorithm.
    716     ///
    717     /// \return <tt>(*this)</tt>
    718     template <typename LowerMap>
    719     NetworkSimplex& lowerMap(const LowerMap& map) {
    720       _have_lower = true;
    721       for (ArcIt a(_graph); a != INVALID; ++a) {
    722         _lower[_arc_id[a]] = map[a];
    723       }
    724       return *this;
    725     }
    726 
    727     /// \brief Set the upper bounds (capacities) on the arcs.
    728     ///
    729     /// This function sets the upper bounds (capacities) on the arcs.
    730     /// If it is not used before calling \ref run(), the upper bounds
    731     /// will be set to \ref INF on all arcs (i.e. the flow value will be
    732     /// unbounded from above on each arc).
    733     ///
    734     /// \param map An arc map storing the upper bounds.
    735     /// Its \c Value type must be convertible to the \c Value type
    736     /// of the algorithm.
    737     ///
    738     /// \return <tt>(*this)</tt>
    739     template<typename UpperMap>
    740     NetworkSimplex& upperMap(const UpperMap& map) {
    741       for (ArcIt a(_graph); a != INVALID; ++a) {
    742         _upper[_arc_id[a]] = map[a];
    743       }
    744       return *this;
    745     }
    746 
    747     /// \brief Set the costs of the arcs.
    748     ///
    749     /// This function sets the costs of the arcs.
    750     /// If it is not used before calling \ref run(), the costs
    751     /// will be set to \c 1 on all arcs.
    752     ///
    753     /// \param map An arc map storing the costs.
    754     /// Its \c Value type must be convertible to the \c Cost type
    755     /// of the algorithm.
    756     ///
    757     /// \return <tt>(*this)</tt>
    758     template<typename CostMap>
    759     NetworkSimplex& costMap(const CostMap& map) {
    760       for (ArcIt a(_graph); a != INVALID; ++a) {
    761         _cost[_arc_id[a]] = map[a];
    762       }
    763       return *this;
    764     }
    765 
    766     /// \brief Set the supply values of the nodes.
    767     ///
    768     /// This function sets the supply values of the nodes.
    769     /// If neither this function nor \ref stSupply() is used before
    770     /// calling \ref run(), the supply of each node will be set to zero.
    771     /// (It makes sense only if non-zero lower bounds are given.)
    772     ///
    773     /// \param map A node map storing the supply values.
    774     /// Its \c Value type must be convertible to the \c Value type
    775     /// of the algorithm.
    776     ///
    777     /// \return <tt>(*this)</tt>
    778     template<typename SupplyMap>
    779     NetworkSimplex& supplyMap(const SupplyMap& map) {
    780       for (NodeIt n(_graph); n != INVALID; ++n) {
    781         _supply[_node_id[n]] = map[n];
    782       }
    783       return *this;
    784     }
    785 
    786     /// \brief Set single source and target nodes and a supply value.
    787     ///
    788     /// This function sets a single source node and a single target node
    789     /// and the required flow value.
    790     /// If neither this function nor \ref supplyMap() is used before
    791     /// calling \ref run(), the supply of each node will be set to zero.
    792     /// (It makes sense only if non-zero lower bounds are given.)
    793     ///
    794     /// Using this function has the same effect as using \ref supplyMap()
    795     /// with such a map in which \c k is assigned to \c s, \c -k is
    796     /// assigned to \c t and all other nodes have zero supply value.
    797     ///
    798     /// \param s The source node.
    799     /// \param t The target node.
    800     /// \param k The required amount of flow from node \c s to node \c t
    801     /// (i.e. the supply of \c s and the demand of \c t).
    802     ///
    803     /// \return <tt>(*this)</tt>
    804     NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
    805       for (int i = 0; i != _node_num; ++i) {
    806         _supply[i] = 0;
    807       }
    808       _supply[_node_id[s]] =  k;
    809       _supply[_node_id[t]] = -k;
    810       return *this;
    811     }
    812    
    813     /// \brief Set the type of the supply constraints.
    814     ///
    815     /// This function sets the type of the supply/demand constraints.
    816     /// If it is not used before calling \ref run(), the \ref GEQ supply
    817     /// type will be used.
    818     ///
    819     /// For more information see \ref SupplyType.
    820     ///
    821     /// \return <tt>(*this)</tt>
    822     NetworkSimplex& supplyType(SupplyType supply_type) {
    823       _stype = supply_type;
    824       return *this;
    825     }
    826 
    827     /// @}
    828 
    829     /// \name Execution Control
    830     /// The algorithm can be executed using \ref run().
    831 
    832     /// @{
    833 
    834     /// \brief Run the algorithm.
    835     ///
    836     /// This function runs the algorithm.
    837     /// The paramters can be specified using functions \ref lowerMap(),
    838     /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
    839     /// \ref supplyType().
    840     /// For example,
    841     /// \code
    842     ///   NetworkSimplex<ListDigraph> ns(graph);
    843     ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
    844     ///     .supplyMap(sup).run();
    845     /// \endcode
    846     ///
    847     /// This function can be called more than once. All the parameters
    848     /// that have been given are kept for the next call, unless
    849     /// \ref reset() is called, thus only the modified parameters
    850     /// have to be set again. See \ref reset() for examples.
    851     /// However the underlying digraph must not be modified after this
    852     /// class have been constructed, since it copies and extends the graph.
    853     ///
    854     /// \param pivot_rule The pivot rule that will be used during the
    855     /// algorithm. For more information see \ref PivotRule.
    856     ///
    857     /// \return \c INFEASIBLE if no feasible flow exists,
    858     /// \n \c OPTIMAL if the problem has optimal solution
    859     /// (i.e. it is feasible and bounded), and the algorithm has found
    860     /// optimal flow and node potentials (primal and dual solutions),
    861     /// \n \c UNBOUNDED if the objective function of the problem is
    862     /// unbounded, i.e. there is a directed cycle having negative total
    863     /// cost and infinite upper bound.
    864     ///
    865     /// \see ProblemType, PivotRule
    866     ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
    867       if (!init()) return INFEASIBLE;
    868       return start(pivot_rule);
    869     }
    870 
    871     /// \brief Reset all the parameters that have been given before.
    872     ///
    873     /// This function resets all the paramaters that have been given
    874     /// before using functions \ref lowerMap(), \ref upperMap(),
    875     /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
    876     ///
    877     /// It is useful for multiple run() calls. If this function is not
    878     /// used, all the parameters given before are kept for the next
    879     /// \ref run() call.
    880     /// However the underlying digraph must not be modified after this
    881     /// class have been constructed, since it copies and extends the graph.
    882     ///
    883     /// For example,
    884     /// \code
    885     ///   NetworkSimplex<ListDigraph> ns(graph);
    886     ///
    887     ///   // First run
    888     ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
    889     ///     .supplyMap(sup).run();
    890     ///
    891     ///   // Run again with modified cost map (reset() is not called,
    892     ///   // so only the cost map have to be set again)
    893     ///   cost[e] += 100;
    894     ///   ns.costMap(cost).run();
    895     ///
    896     ///   // Run again from scratch using reset()
    897     ///   // (the lower bounds will be set to zero on all arcs)
    898     ///   ns.reset();
    899     ///   ns.upperMap(capacity).costMap(cost)
    900     ///     .supplyMap(sup).run();
    901     /// \endcode
    902     ///
    903     /// \return <tt>(*this)</tt>
    904     NetworkSimplex& reset() {
    905       for (int i = 0; i != _node_num; ++i) {
    906         _supply[i] = 0;
    907       }
    908       for (int i = 0; i != _arc_num; ++i) {
    909         _lower[i] = 0;
    910         _upper[i] = INF;
    911         _cost[i] = 1;
    912       }
    913       _have_lower = false;
    914       _stype = GEQ;
     928      if (_arc_mixing) {
     929        // Store the arcs in a mixed order
     930        int k = std::max(int(std::sqrt(double(_arc_num))), 10);
     931        int i = 0, j = 0;
     932        for (ArcIt a(_graph); a != INVALID; ++a) {
     933          _arc_id[a] = i;
     934          _source[i] = _node_id[_graph.source(a)];
     935          _target[i] = _node_id[_graph.target(a)];
     936          if ((i += k) >= _arc_num) i = ++j;
     937        }
     938      } else {
     939        // Store the arcs in the original order
     940        int i = 0;
     941        for (ArcIt a(_graph); a != INVALID; ++a, ++i) {
     942          _arc_id[a] = i;
     943          _source[i] = _node_id[_graph.source(a)];
     944          _target[i] = _node_id[_graph.target(a)];
     945        }
     946      }
     947
     948      // Reset parameters
     949      resetParams();
    915950      return *this;
    916951    }
     
    10251060          Value c = _lower[i];
    10261061          if (c >= 0) {
    1027             _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
     1062            _cap[i] = _upper[i] < MAX ? _upper[i] - c : INF;
    10281063          } else {
    1029             _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
     1064            _cap[i] = _upper[i] < MAX + c ? _upper[i] - c : INF;
    10301065          }
    10311066          _supply[_source[i]] -= c;
     
    10551090        _state[i] = STATE_LOWER;
    10561091      }
    1057      
     1092
    10581093      // Set data for the artificial root node
    10591094      _root = _node_num;
     
    12191254        e = _pred[u];
    12201255        d = _forward[u] ?
    1221           _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
     1256          _flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]);
    12221257        if (d < delta) {
    12231258          delta = d;
     
    12291264      for (int u = second; u != join; u = _parent[u]) {
    12301265        e = _pred[u];
    1231         d = _forward[u] ? 
    1232           (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
     1266        d = _forward[u] ?
     1267          (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
    12331268        if (d <= delta) {
    12341269          delta = d;
     
    13311366
    13321367      // Update _rev_thread using the new _thread values
    1333       for (int i = 0; i < int(_dirty_revs.size()); ++i) {
     1368      for (int i = 0; i != int(_dirty_revs.size()); ++i) {
    13341369        u = _dirty_revs[i];
    13351370        _rev_thread[_thread[u]] = u;
     
    14011436        _pi[u] += sigma;
    14021437      }
     1438    }
     1439
     1440    // Heuristic initial pivots
     1441    bool initialPivots() {
     1442      Value curr, total = 0;
     1443      std::vector<Node> supply_nodes, demand_nodes;
     1444      for (NodeIt u(_graph); u != INVALID; ++u) {
     1445        curr = _supply[_node_id[u]];
     1446        if (curr > 0) {
     1447          total += curr;
     1448          supply_nodes.push_back(u);
     1449        }
     1450        else if (curr < 0) {
     1451          demand_nodes.push_back(u);
     1452        }
     1453      }
     1454      if (_sum_supply > 0) total -= _sum_supply;
     1455      if (total <= 0) return true;
     1456
     1457      IntVector arc_vector;
     1458      if (_sum_supply >= 0) {
     1459        if (supply_nodes.size() == 1 && demand_nodes.size() == 1) {
     1460          // Perform a reverse graph search from the sink to the source
     1461          typename GR::template NodeMap<bool> reached(_graph, false);
     1462          Node s = supply_nodes[0], t = demand_nodes[0];
     1463          std::vector<Node> stack;
     1464          reached[t] = true;
     1465          stack.push_back(t);
     1466          while (!stack.empty()) {
     1467            Node u, v = stack.back();
     1468            stack.pop_back();
     1469            if (v == s) break;
     1470            for (InArcIt a(_graph, v); a != INVALID; ++a) {
     1471              if (reached[u = _graph.source(a)]) continue;
     1472              int j = _arc_id[a];
     1473              if (_cap[j] >= total) {
     1474                arc_vector.push_back(j);
     1475                reached[u] = true;
     1476                stack.push_back(u);
     1477              }
     1478            }
     1479          }
     1480        } else {
     1481          // Find the min. cost incomming arc for each demand node
     1482          for (int i = 0; i != int(demand_nodes.size()); ++i) {
     1483            Node v = demand_nodes[i];
     1484            Cost c, min_cost = std::numeric_limits<Cost>::max();
     1485            Arc min_arc = INVALID;
     1486            for (InArcIt a(_graph, v); a != INVALID; ++a) {
     1487              c = _cost[_arc_id[a]];
     1488              if (c < min_cost) {
     1489                min_cost = c;
     1490                min_arc = a;
     1491              }
     1492            }
     1493            if (min_arc != INVALID) {
     1494              arc_vector.push_back(_arc_id[min_arc]);
     1495            }
     1496          }
     1497        }
     1498      } else {
     1499        // Find the min. cost outgoing arc for each supply node
     1500        for (int i = 0; i != int(supply_nodes.size()); ++i) {
     1501          Node u = supply_nodes[i];
     1502          Cost c, min_cost = std::numeric_limits<Cost>::max();
     1503          Arc min_arc = INVALID;
     1504          for (OutArcIt a(_graph, u); a != INVALID; ++a) {
     1505            c = _cost[_arc_id[a]];
     1506            if (c < min_cost) {
     1507              min_cost = c;
     1508              min_arc = a;
     1509            }
     1510          }
     1511          if (min_arc != INVALID) {
     1512            arc_vector.push_back(_arc_id[min_arc]);
     1513          }
     1514        }
     1515      }
     1516
     1517      // Perform heuristic initial pivots
     1518      for (int i = 0; i != int(arc_vector.size()); ++i) {
     1519        in_arc = arc_vector[i];
     1520        if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] -
     1521            _pi[_target[in_arc]]) >= 0) continue;
     1522        findJoinNode();
     1523        bool change = findLeavingArc();
     1524        if (delta >= MAX) return false;
     1525        changeFlow(change);
     1526        if (change) {
     1527          updateTreeStructure();
     1528          updatePotential();
     1529        }
     1530      }
     1531      return true;
    14031532    }
    14041533
     
    14251554      PivotRuleImpl pivot(*this);
    14261555
     1556      // Perform heuristic initial pivots
     1557      if (!initialPivots()) return UNBOUNDED;
     1558
    14271559      // Execute the Network Simplex algorithm
    14281560      while (pivot.findEnteringArc()) {
    14291561        findJoinNode();
    14301562        bool change = findLeavingArc();
    1431         if (delta >= INF) return UNBOUNDED;
     1563        if (delta >= MAX) return UNBOUNDED;
    14321564        changeFlow(change);
    14331565        if (change) {
     
    14361568        }
    14371569      }
    1438      
     1570
    14391571      // Check feasibility
    14401572      for (int e = _search_arc_num; e != _all_arc_num; ++e) {
     
    14531585        }
    14541586      }
    1455      
     1587
    14561588      // Shift potentials to meet the requirements of the GEQ/LEQ type
    14571589      // optimality conditions
Note: See TracChangeset for help on using the changeset viewer.