COIN-OR::LEMON - Graph Library

source: lemon/lemon/network_simplex.h @ 1270:dceba191c00d

Last change on this file since 1270:dceba191c00d was 1270:dceba191c00d, checked in by Alpar Juttner <alpar@…>, 11 years ago

Apply unify-sources.sh to the source tree

File size: 51.6 KB
RevLine 
[648]1/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library.
4 *
[1270]5 * Copyright (C) 2003-2013
[648]6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
12 *
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_NETWORK_SIMPLEX_H
20#define LEMON_NETWORK_SIMPLEX_H
21
[710]22/// \ingroup min_cost_flow_algs
[648]23///
24/// \file
[652]25/// \brief Network Simplex algorithm for finding a minimum cost flow.
[648]26
27#include <vector>
28#include <limits>
29#include <algorithm>
30
[650]31#include <lemon/core.h>
[648]32#include <lemon/math.h>
33
34namespace lemon {
35
[710]36  /// \addtogroup min_cost_flow_algs
[648]37  /// @{
38
[652]39  /// \brief Implementation of the primal Network Simplex algorithm
[648]40  /// for finding a \ref min_cost_flow "minimum cost flow".
41  ///
[652]42  /// \ref NetworkSimplex implements the primal Network Simplex algorithm
[802]43  /// for finding a \ref min_cost_flow "minimum cost flow"
[1221]44  /// \cite amo93networkflows, \cite dantzig63linearprog,
45  /// \cite kellyoneill91netsimplex.
[878]46  /// This algorithm is a highly efficient specialized version of the
47  /// linear programming simplex method directly for the minimum cost
48  /// flow problem.
[653]49  ///
[1023]50  /// In general, \ref NetworkSimplex and \ref CostScaling are the fastest
[1165]51  /// implementations available in LEMON for solving this problem.
52  /// (For more information, see \ref min_cost_flow_algs "the module page".)
[1023]53  /// Furthermore, this class supports both directions of the supply/demand
54  /// inequality constraints. For more information, see \ref SupplyType.
[687]55  ///
56  /// Most of the parameters of the problem (except for the digraph)
57  /// can be given using separate functions, and the algorithm can be
58  /// executed using the \ref run() function. If some parameters are not
59  /// specified, then default values will be used.
[648]60  ///
[652]61  /// \tparam GR The digraph type the algorithm runs on.
[878]62  /// \tparam V The number type used for flow amounts, capacity bounds
[833]63  /// and supply values in the algorithm. By default, it is \c int.
[878]64  /// \tparam C The number type used for costs and potentials in the
[833]65  /// algorithm. By default, it is the same as \c V.
[648]66  ///
[1025]67  /// \warning Both \c V and \c C must be signed number types.
68  /// \warning All input data (capacities, supply values, and costs) must
[655]69  /// be integer.
[648]70  ///
[652]71  /// \note %NetworkSimplex provides five different pivot rule
[656]72  /// implementations, from which the most efficient one is used
[833]73  /// by default. For more information, see \ref PivotRule.
[688]74  template <typename GR, typename V = int, typename C = V>
[648]75  class NetworkSimplex
76  {
[652]77  public:
[648]78
[689]79    /// The type of the flow amounts, capacity bounds and supply values
[688]80    typedef V Value;
[689]81    /// The type of the arc costs
[654]82    typedef C Cost;
[652]83
84  public:
85
[687]86    /// \brief Problem type constants for the \c run() function.
[652]87    ///
[687]88    /// Enum type containing the problem type constants that can be
89    /// returned by the \ref run() function of the algorithm.
90    enum ProblemType {
91      /// The problem has no feasible solution (flow).
92      INFEASIBLE,
93      /// The problem has optimal solution (i.e. it is feasible and
94      /// bounded), and the algorithm has found optimal flow and node
95      /// potentials (primal and dual solutions).
96      OPTIMAL,
97      /// The objective function of the problem is unbounded, i.e.
98      /// there is a directed cycle having negative total cost and
99      /// infinite upper bound.
100      UNBOUNDED
101    };
[956]102
[687]103    /// \brief Constants for selecting the type of the supply constraints.
104    ///
105    /// Enum type containing constants for selecting the supply type,
106    /// i.e. the direction of the inequalities in the supply/demand
107    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
108    ///
[710]109    /// The default supply type is \c GEQ, the \c LEQ type can be
110    /// selected using \ref supplyType().
111    /// The equality form is a special case of both supply types.
[687]112    enum SupplyType {
113      /// This option means that there are <em>"greater or equal"</em>
[710]114      /// supply/demand constraints in the definition of the problem.
[687]115      GEQ,
116      /// This option means that there are <em>"less or equal"</em>
[710]117      /// supply/demand constraints in the definition of the problem.
118      LEQ
[687]119    };
[956]120
[687]121    /// \brief Constants for selecting the pivot rule.
122    ///
123    /// Enum type containing constants for selecting the pivot rule for
124    /// the \ref run() function.
125    ///
[1136]126    /// \ref NetworkSimplex provides five different implementations for
127    /// the pivot strategy that significantly affects the running time
[652]128    /// of the algorithm.
[1136]129    /// According to experimental tests conducted on various problem
130    /// instances, \ref BLOCK_SEARCH "Block Search" and
131    /// \ref ALTERING_LIST "Altering Candidate List" rules turned out
132    /// to be the most efficient.
133    /// Since \ref BLOCK_SEARCH "Block Search" is a simpler strategy that
134    /// seemed to be slightly more robust, it is used by default.
135    /// However, another pivot rule can easily be selected using the
136    /// \ref run() function with the proper parameter.
[652]137    enum PivotRule {
138
[833]139      /// The \e First \e Eligible pivot rule.
[652]140      /// The next eligible arc is selected in a wraparound fashion
141      /// in every iteration.
142      FIRST_ELIGIBLE,
143
[833]144      /// The \e Best \e Eligible pivot rule.
[652]145      /// The best eligible arc is selected in every iteration.
146      BEST_ELIGIBLE,
147
[833]148      /// The \e Block \e Search pivot rule.
[652]149      /// A specified number of arcs are examined in every iteration
150      /// in a wraparound fashion and the best eligible arc is selected
151      /// from this block.
152      BLOCK_SEARCH,
153
[833]154      /// The \e Candidate \e List pivot rule.
[652]155      /// In a major iteration a candidate list is built from eligible arcs
156      /// in a wraparound fashion and in the following minor iterations
157      /// the best eligible arc is selected from this list.
158      CANDIDATE_LIST,
159
[833]160      /// The \e Altering \e Candidate \e List pivot rule.
[652]161      /// It is a modified version of the Candidate List method.
[1136]162      /// It keeps only a few of the best eligible arcs from the former
[652]163      /// candidate list and extends this list in every iteration.
164      ALTERING_LIST
165    };
[956]166
[652]167  private:
168
169    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
170
[648]171    typedef std::vector<int> IntVector;
[689]172    typedef std::vector<Value> ValueVector;
[654]173    typedef std::vector<Cost> CostVector;
[990]174    typedef std::vector<signed char> CharVector;
[1023]175    // Note: vector<signed char> is used instead of vector<ArcState> and
[990]176    // vector<ArcDirection> for efficiency reasons
[648]177
178    // State constants for arcs
[936]179    enum ArcState {
[648]180      STATE_UPPER = -1,
181      STATE_TREE  =  0,
182      STATE_LOWER =  1
183    };
184
[990]185    // Direction constants for tree arcs
186    enum ArcDirection {
187      DIR_DOWN = -1,
188      DIR_UP   =  1
189    };
[936]190
[648]191  private:
192
[652]193    // Data related to the underlying digraph
194    const GR &_graph;
195    int _node_num;
196    int _arc_num;
[710]197    int _all_arc_num;
198    int _search_arc_num;
[652]199
200    // Parameters of the problem
[689]201    bool _have_lower;
[687]202    SupplyType _stype;
[688]203    Value _sum_supply;
[648]204
[652]205    // Data structures for storing the digraph
[650]206    IntNodeMap _node_id;
[689]207    IntArcMap _arc_id;
[650]208    IntVector _source;
209    IntVector _target;
[898]210    bool _arc_mixing;
[650]211
[652]212    // Node and arc data
[689]213    ValueVector _lower;
214    ValueVector _upper;
215    ValueVector _cap;
[654]216    CostVector _cost;
[689]217    ValueVector _supply;
218    ValueVector _flow;
[654]219    CostVector _pi;
[648]220
[650]221    // Data for storing the spanning tree structure
[648]222    IntVector _parent;
223    IntVector _pred;
224    IntVector _thread;
[651]225    IntVector _rev_thread;
226    IntVector _succ_num;
227    IntVector _last_succ;
[990]228    CharVector _pred_dir;
229    CharVector _state;
[651]230    IntVector _dirty_revs;
[648]231    int _root;
232
233    // Temporary data used in the current pivot iteration
[650]234    int in_arc, join, u_in, v_in, u_out, v_out;
[688]235    Value delta;
[648]236
[877]237    const Value MAX;
[710]238
[687]239  public:
[956]240
[687]241    /// \brief Constant for infinite upper bounds (capacities).
242    ///
243    /// Constant for infinite upper bounds (capacities).
[688]244    /// It is \c std::numeric_limits<Value>::infinity() if available,
245    /// \c std::numeric_limits<Value>::max() otherwise.
246    const Value INF;
[687]247
[648]248  private:
249
[652]250    // Implementation of the First Eligible pivot rule
[648]251    class FirstEligiblePivotRule
252    {
253    private:
254
255      // References to the NetworkSimplex class
256      const IntVector  &_source;
257      const IntVector  &_target;
[654]258      const CostVector &_cost;
[990]259      const CharVector &_state;
[654]260      const CostVector &_pi;
[648]261      int &_in_arc;
[710]262      int _search_arc_num;
[648]263
264      // Pivot rule data
265      int _next_arc;
266
267    public:
268
[652]269      // Constructor
[648]270      FirstEligiblePivotRule(NetworkSimplex &ns) :
[650]271        _source(ns._source), _target(ns._target),
[648]272        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
[710]273        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
274        _next_arc(0)
[648]275      {}
276
[652]277      // Find next entering arc
[648]278      bool findEnteringArc() {
[654]279        Cost c;
[910]280        for (int e = _next_arc; e != _search_arc_num; ++e) {
[648]281          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
282          if (c < 0) {
283            _in_arc = e;
284            _next_arc = e + 1;
285            return true;
286          }
287        }
[910]288        for (int e = 0; e != _next_arc; ++e) {
[648]289          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
290          if (c < 0) {
291            _in_arc = e;
292            _next_arc = e + 1;
293            return true;
294          }
295        }
296        return false;
297      }
298
299    }; //class FirstEligiblePivotRule
300
301
[652]302    // Implementation of the Best Eligible pivot rule
[648]303    class BestEligiblePivotRule
304    {
305    private:
306
307      // References to the NetworkSimplex class
308      const IntVector  &_source;
309      const IntVector  &_target;
[654]310      const CostVector &_cost;
[990]311      const CharVector &_state;
[654]312      const CostVector &_pi;
[648]313      int &_in_arc;
[710]314      int _search_arc_num;
[648]315
316    public:
317
[652]318      // Constructor
[648]319      BestEligiblePivotRule(NetworkSimplex &ns) :
[650]320        _source(ns._source), _target(ns._target),
[648]321        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
[710]322        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num)
[648]323      {}
324
[652]325      // Find next entering arc
[648]326      bool findEnteringArc() {
[654]327        Cost c, min = 0;
[910]328        for (int e = 0; e != _search_arc_num; ++e) {
[648]329          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
330          if (c < min) {
331            min = c;
332            _in_arc = e;
333          }
334        }
335        return min < 0;
336      }
337
338    }; //class BestEligiblePivotRule
339
340
[652]341    // Implementation of the Block Search pivot rule
[648]342    class BlockSearchPivotRule
343    {
344    private:
345
346      // References to the NetworkSimplex class
347      const IntVector  &_source;
348      const IntVector  &_target;
[654]349      const CostVector &_cost;
[990]350      const CharVector &_state;
[654]351      const CostVector &_pi;
[648]352      int &_in_arc;
[710]353      int _search_arc_num;
[648]354
355      // Pivot rule data
356      int _block_size;
357      int _next_arc;
358
359    public:
360
[652]361      // Constructor
[648]362      BlockSearchPivotRule(NetworkSimplex &ns) :
[650]363        _source(ns._source), _target(ns._target),
[648]364        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
[710]365        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
366        _next_arc(0)
[648]367      {
368        // The main parameters of the pivot rule
[910]369        const double BLOCK_SIZE_FACTOR = 1.0;
[648]370        const int MIN_BLOCK_SIZE = 10;
371
[659]372        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
[710]373                                    std::sqrt(double(_search_arc_num))),
[648]374                                MIN_BLOCK_SIZE );
375      }
376
[652]377      // Find next entering arc
[648]378      bool findEnteringArc() {
[654]379        Cost c, min = 0;
[648]380        int cnt = _block_size;
[774]381        int e;
[910]382        for (e = _next_arc; e != _search_arc_num; ++e) {
[648]383          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
384          if (c < min) {
385            min = c;
[774]386            _in_arc = e;
[648]387          }
388          if (--cnt == 0) {
[774]389            if (min < 0) goto search_end;
[648]390            cnt = _block_size;
391          }
392        }
[910]393        for (e = 0; e != _next_arc; ++e) {
[774]394          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
395          if (c < min) {
396            min = c;
397            _in_arc = e;
398          }
399          if (--cnt == 0) {
400            if (min < 0) goto search_end;
401            cnt = _block_size;
[648]402          }
403        }
404        if (min >= 0) return false;
[774]405
406      search_end:
[648]407        _next_arc = e;
408        return true;
409      }
410
411    }; //class BlockSearchPivotRule
412
413
[652]414    // Implementation of the Candidate List pivot rule
[648]415    class CandidateListPivotRule
416    {
417    private:
418
419      // References to the NetworkSimplex class
420      const IntVector  &_source;
421      const IntVector  &_target;
[654]422      const CostVector &_cost;
[990]423      const CharVector &_state;
[654]424      const CostVector &_pi;
[648]425      int &_in_arc;
[710]426      int _search_arc_num;
[648]427
428      // Pivot rule data
429      IntVector _candidates;
430      int _list_length, _minor_limit;
431      int _curr_length, _minor_count;
432      int _next_arc;
433
434    public:
435
436      /// Constructor
437      CandidateListPivotRule(NetworkSimplex &ns) :
[650]438        _source(ns._source), _target(ns._target),
[648]439        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
[710]440        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
441        _next_arc(0)
[648]442      {
443        // The main parameters of the pivot rule
[774]444        const double LIST_LENGTH_FACTOR = 0.25;
[648]445        const int MIN_LIST_LENGTH = 10;
446        const double MINOR_LIMIT_FACTOR = 0.1;
447        const int MIN_MINOR_LIMIT = 3;
448
[659]449        _list_length = std::max( int(LIST_LENGTH_FACTOR *
[710]450                                     std::sqrt(double(_search_arc_num))),
[648]451                                 MIN_LIST_LENGTH );
452        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
453                                 MIN_MINOR_LIMIT );
454        _curr_length = _minor_count = 0;
455        _candidates.resize(_list_length);
456      }
457
458      /// Find next entering arc
459      bool findEnteringArc() {
[654]460        Cost min, c;
[774]461        int e;
[648]462        if (_curr_length > 0 && _minor_count < _minor_limit) {
463          // Minor iteration: select the best eligible arc from the
464          // current candidate list
465          ++_minor_count;
466          min = 0;
467          for (int i = 0; i < _curr_length; ++i) {
468            e = _candidates[i];
469            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
470            if (c < min) {
471              min = c;
[774]472              _in_arc = e;
[648]473            }
[774]474            else if (c >= 0) {
[648]475              _candidates[i--] = _candidates[--_curr_length];
476            }
477          }
[774]478          if (min < 0) return true;
[648]479        }
480
481        // Major iteration: build a new candidate list
482        min = 0;
483        _curr_length = 0;
[910]484        for (e = _next_arc; e != _search_arc_num; ++e) {
[648]485          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
486          if (c < 0) {
487            _candidates[_curr_length++] = e;
488            if (c < min) {
489              min = c;
[774]490              _in_arc = e;
[648]491            }
[774]492            if (_curr_length == _list_length) goto search_end;
[648]493          }
494        }
[910]495        for (e = 0; e != _next_arc; ++e) {
[774]496          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
497          if (c < 0) {
498            _candidates[_curr_length++] = e;
499            if (c < min) {
500              min = c;
501              _in_arc = e;
[648]502            }
[774]503            if (_curr_length == _list_length) goto search_end;
[648]504          }
505        }
506        if (_curr_length == 0) return false;
[956]507
508      search_end:
[648]509        _minor_count = 1;
510        _next_arc = e;
511        return true;
512      }
513
514    }; //class CandidateListPivotRule
515
516
[652]517    // Implementation of the Altering Candidate List pivot rule
[648]518    class AlteringListPivotRule
519    {
520    private:
521
522      // References to the NetworkSimplex class
523      const IntVector  &_source;
524      const IntVector  &_target;
[654]525      const CostVector &_cost;
[990]526      const CharVector &_state;
[654]527      const CostVector &_pi;
[648]528      int &_in_arc;
[710]529      int _search_arc_num;
[648]530
531      // Pivot rule data
532      int _block_size, _head_length, _curr_length;
533      int _next_arc;
534      IntVector _candidates;
[654]535      CostVector _cand_cost;
[648]536
537      // Functor class to compare arcs during sort of the candidate list
538      class SortFunc
539      {
540      private:
[654]541        const CostVector &_map;
[648]542      public:
[654]543        SortFunc(const CostVector &map) : _map(map) {}
[648]544        bool operator()(int left, int right) {
[1136]545          return _map[left] < _map[right];
[648]546        }
547      };
548
549      SortFunc _sort_func;
550
551    public:
552
[652]553      // Constructor
[648]554      AlteringListPivotRule(NetworkSimplex &ns) :
[650]555        _source(ns._source), _target(ns._target),
[648]556        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
[710]557        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
558        _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
[648]559      {
560        // The main parameters of the pivot rule
[774]561        const double BLOCK_SIZE_FACTOR = 1.0;
[648]562        const int MIN_BLOCK_SIZE = 10;
[1136]563        const double HEAD_LENGTH_FACTOR = 0.01;
[648]564        const int MIN_HEAD_LENGTH = 3;
565
[659]566        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
[710]567                                    std::sqrt(double(_search_arc_num))),
[648]568                                MIN_BLOCK_SIZE );
569        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
570                                 MIN_HEAD_LENGTH );
571        _candidates.resize(_head_length + _block_size);
572        _curr_length = 0;
573      }
574
[652]575      // Find next entering arc
[648]576      bool findEnteringArc() {
577        // Check the current candidate list
578        int e;
[990]579        Cost c;
[910]580        for (int i = 0; i != _curr_length; ++i) {
[648]581          e = _candidates[i];
[990]582          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
583          if (c < 0) {
584            _cand_cost[e] = c;
585          } else {
[648]586            _candidates[i--] = _candidates[--_curr_length];
587          }
588        }
589
590        // Extend the list
591        int cnt = _block_size;
592        int limit = _head_length;
593
[910]594        for (e = _next_arc; e != _search_arc_num; ++e) {
[990]595          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
596          if (c < 0) {
597            _cand_cost[e] = c;
[648]598            _candidates[_curr_length++] = e;
599          }
600          if (--cnt == 0) {
[774]601            if (_curr_length > limit) goto search_end;
[648]602            limit = 0;
603            cnt = _block_size;
604          }
605        }
[910]606        for (e = 0; e != _next_arc; ++e) {
[1136]607          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
608          if (c < 0) {
609            _cand_cost[e] = c;
[774]610            _candidates[_curr_length++] = e;
611          }
612          if (--cnt == 0) {
613            if (_curr_length > limit) goto search_end;
614            limit = 0;
615            cnt = _block_size;
[648]616          }
617        }
618        if (_curr_length == 0) return false;
[956]619
[774]620      search_end:
[648]621
[1136]622        // Perform partial sort operation on the candidate list
623        int new_length = std::min(_head_length + 1, _curr_length);
624        std::partial_sort(_candidates.begin(), _candidates.begin() + new_length,
625                          _candidates.begin() + _curr_length, _sort_func);
[648]626
[1136]627        // Select the entering arc and remove it from the list
[648]628        _in_arc = _candidates[0];
[774]629        _next_arc = e;
[1136]630        _candidates[0] = _candidates[new_length - 1];
631        _curr_length = new_length - 1;
[648]632        return true;
633      }
634
635    }; //class AlteringListPivotRule
636
637  public:
638
[652]639    /// \brief Constructor.
[648]640    ///
[656]641    /// The constructor of the class.
[648]642    ///
[650]643    /// \param graph The digraph the algorithm runs on.
[991]644    /// \param arc_mixing Indicate if the arcs will be stored in a
[956]645    /// mixed order in the internal data structure.
[991]646    /// In general, it leads to similar performance as using the original
647    /// arc order, but it makes the algorithm more robust and in special
648    /// cases, even significantly faster. Therefore, it is enabled by default.
649    NetworkSimplex(const GR& graph, bool arc_mixing = true) :
[689]650      _graph(graph), _node_id(graph), _arc_id(graph),
[898]651      _arc_mixing(arc_mixing),
[877]652      MAX(std::numeric_limits<Value>::max()),
[688]653      INF(std::numeric_limits<Value>::has_infinity ?
[877]654          std::numeric_limits<Value>::infinity() : MAX)
[652]655    {
[878]656      // Check the number types
[688]657      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
[687]658        "The flow type of NetworkSimplex must be signed");
659      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
660        "The cost type of NetworkSimplex must be signed");
[648]661
[898]662      // Reset data structures
[776]663      reset();
[648]664    }
665
[656]666    /// \name Parameters
667    /// The parameters of the algorithm can be specified using these
668    /// functions.
669
670    /// @{
671
[652]672    /// \brief Set the lower bounds on the arcs.
673    ///
674    /// This function sets the lower bounds on the arcs.
[687]675    /// If it is not used before calling \ref run(), the lower bounds
676    /// will be set to zero on all arcs.
[652]677    ///
678    /// \param map An arc map storing the lower bounds.
[688]679    /// Its \c Value type must be convertible to the \c Value type
[652]680    /// of the algorithm.
681    ///
682    /// \return <tt>(*this)</tt>
[687]683    template <typename LowerMap>
684    NetworkSimplex& lowerMap(const LowerMap& map) {
[689]685      _have_lower = true;
[652]686      for (ArcIt a(_graph); a != INVALID; ++a) {
[689]687        _lower[_arc_id[a]] = map[a];
[652]688      }
689      return *this;
690    }
691
692    /// \brief Set the upper bounds (capacities) on the arcs.
693    ///
694    /// This function sets the upper bounds (capacities) on the arcs.
[687]695    /// If it is not used before calling \ref run(), the upper bounds
696    /// will be set to \ref INF on all arcs (i.e. the flow value will be
[878]697    /// unbounded from above).
[652]698    ///
699    /// \param map An arc map storing the upper bounds.
[688]700    /// Its \c Value type must be convertible to the \c Value type
[652]701    /// of the algorithm.
702    ///
703    /// \return <tt>(*this)</tt>
[687]704    template<typename UpperMap>
705    NetworkSimplex& upperMap(const UpperMap& map) {
[652]706      for (ArcIt a(_graph); a != INVALID; ++a) {
[689]707        _upper[_arc_id[a]] = map[a];
[652]708      }
709      return *this;
710    }
711
712    /// \brief Set the costs of the arcs.
713    ///
714    /// This function sets the costs of the arcs.
715    /// If it is not used before calling \ref run(), the costs
716    /// will be set to \c 1 on all arcs.
717    ///
718    /// \param map An arc map storing the costs.
[654]719    /// Its \c Value type must be convertible to the \c Cost type
[652]720    /// of the algorithm.
721    ///
722    /// \return <tt>(*this)</tt>
[687]723    template<typename CostMap>
724    NetworkSimplex& costMap(const CostMap& map) {
[652]725      for (ArcIt a(_graph); a != INVALID; ++a) {
[689]726        _cost[_arc_id[a]] = map[a];
[652]727      }
728      return *this;
729    }
730
731    /// \brief Set the supply values of the nodes.
732    ///
733    /// This function sets the supply values of the nodes.
734    /// If neither this function nor \ref stSupply() is used before
735    /// calling \ref run(), the supply of each node will be set to zero.
736    ///
737    /// \param map A node map storing the supply values.
[688]738    /// Its \c Value type must be convertible to the \c Value type
[652]739    /// of the algorithm.
740    ///
741    /// \return <tt>(*this)</tt>
[1023]742    ///
743    /// \sa supplyType()
[687]744    template<typename SupplyMap>
745    NetworkSimplex& supplyMap(const SupplyMap& map) {
[652]746      for (NodeIt n(_graph); n != INVALID; ++n) {
[689]747        _supply[_node_id[n]] = map[n];
[652]748      }
749      return *this;
750    }
751
752    /// \brief Set single source and target nodes and a supply value.
753    ///
754    /// This function sets a single source node and a single target node
755    /// and the required flow value.
756    /// If neither this function nor \ref supplyMap() is used before
757    /// calling \ref run(), the supply of each node will be set to zero.
758    ///
[687]759    /// Using this function has the same effect as using \ref supplyMap()
[1023]760    /// with a map in which \c k is assigned to \c s, \c -k is
[687]761    /// assigned to \c t and all other nodes have zero supply value.
762    ///
[652]763    /// \param s The source node.
764    /// \param t The target node.
765    /// \param k The required amount of flow from node \c s to node \c t
766    /// (i.e. the supply of \c s and the demand of \c t).
767    ///
768    /// \return <tt>(*this)</tt>
[688]769    NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
[689]770      for (int i = 0; i != _node_num; ++i) {
771        _supply[i] = 0;
772      }
773      _supply[_node_id[s]] =  k;
774      _supply[_node_id[t]] = -k;
[652]775      return *this;
776    }
[956]777
[687]778    /// \brief Set the type of the supply constraints.
[656]779    ///
[687]780    /// This function sets the type of the supply/demand constraints.
781    /// If it is not used before calling \ref run(), the \ref GEQ supply
[656]782    /// type will be used.
783    ///
[833]784    /// For more information, see \ref SupplyType.
[656]785    ///
786    /// \return <tt>(*this)</tt>
[687]787    NetworkSimplex& supplyType(SupplyType supply_type) {
788      _stype = supply_type;
[656]789      return *this;
790    }
[652]791
[656]792    /// @}
[648]793
[652]794    /// \name Execution Control
795    /// The algorithm can be executed using \ref run().
796
[648]797    /// @{
798
799    /// \brief Run the algorithm.
800    ///
801    /// This function runs the algorithm.
[656]802    /// The paramters can be specified using functions \ref lowerMap(),
[956]803    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
[689]804    /// \ref supplyType().
[656]805    /// For example,
[652]806    /// \code
807    ///   NetworkSimplex<ListDigraph> ns(graph);
[687]808    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
[652]809    ///     .supplyMap(sup).run();
810    /// \endcode
[648]811    ///
[898]812    /// This function can be called more than once. All the given parameters
813    /// are kept for the next call, unless \ref resetParams() or \ref reset()
814    /// is used, thus only the modified parameters have to be set again.
815    /// If the underlying digraph was also modified after the construction
816    /// of the class (or the last \ref reset() call), then the \ref reset()
817    /// function must be called.
[653]818    ///
[652]819    /// \param pivot_rule The pivot rule that will be used during the
[833]820    /// algorithm. For more information, see \ref PivotRule.
[648]821    ///
[687]822    /// \return \c INFEASIBLE if no feasible flow exists,
823    /// \n \c OPTIMAL if the problem has optimal solution
824    /// (i.e. it is feasible and bounded), and the algorithm has found
825    /// optimal flow and node potentials (primal and dual solutions),
826    /// \n \c UNBOUNDED if the objective function of the problem is
827    /// unbounded, i.e. there is a directed cycle having negative total
828    /// cost and infinite upper bound.
829    ///
830    /// \see ProblemType, PivotRule
[898]831    /// \see resetParams(), reset()
[687]832    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
833      if (!init()) return INFEASIBLE;
834      return start(pivot_rule);
[648]835    }
836
[653]837    /// \brief Reset all the parameters that have been given before.
838    ///
839    /// This function resets all the paramaters that have been given
[656]840    /// before using functions \ref lowerMap(), \ref upperMap(),
[689]841    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
[653]842    ///
[898]843    /// It is useful for multiple \ref run() calls. Basically, all the given
844    /// parameters are kept for the next \ref run() call, unless
845    /// \ref resetParams() or \ref reset() is used.
846    /// If the underlying digraph was also modified after the construction
847    /// of the class or the last \ref reset() call, then the \ref reset()
848    /// function must be used, otherwise \ref resetParams() is sufficient.
[653]849    ///
850    /// For example,
851    /// \code
852    ///   NetworkSimplex<ListDigraph> ns(graph);
853    ///
854    ///   // First run
[687]855    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
[653]856    ///     .supplyMap(sup).run();
857    ///
[898]858    ///   // Run again with modified cost map (resetParams() is not called,
[653]859    ///   // so only the cost map have to be set again)
860    ///   cost[e] += 100;
861    ///   ns.costMap(cost).run();
862    ///
[898]863    ///   // Run again from scratch using resetParams()
[653]864    ///   // (the lower bounds will be set to zero on all arcs)
[898]865    ///   ns.resetParams();
[687]866    ///   ns.upperMap(capacity).costMap(cost)
[653]867    ///     .supplyMap(sup).run();
868    /// \endcode
869    ///
870    /// \return <tt>(*this)</tt>
[898]871    ///
872    /// \see reset(), run()
873    NetworkSimplex& resetParams() {
[689]874      for (int i = 0; i != _node_num; ++i) {
875        _supply[i] = 0;
876      }
877      for (int i = 0; i != _arc_num; ++i) {
878        _lower[i] = 0;
879        _upper[i] = INF;
880        _cost[i] = 1;
881      }
882      _have_lower = false;
[687]883      _stype = GEQ;
[653]884      return *this;
885    }
886
[898]887    /// \brief Reset the internal data structures and all the parameters
888    /// that have been given before.
889    ///
890    /// This function resets the internal data structures and all the
891    /// paramaters that have been given before using functions \ref lowerMap(),
892    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
893    /// \ref supplyType().
894    ///
895    /// It is useful for multiple \ref run() calls. Basically, all the given
896    /// parameters are kept for the next \ref run() call, unless
897    /// \ref resetParams() or \ref reset() is used.
898    /// If the underlying digraph was also modified after the construction
899    /// of the class or the last \ref reset() call, then the \ref reset()
900    /// function must be used, otherwise \ref resetParams() is sufficient.
901    ///
902    /// See \ref resetParams() for examples.
903    ///
904    /// \return <tt>(*this)</tt>
905    ///
906    /// \see resetParams(), run()
907    NetworkSimplex& reset() {
908      // Resize vectors
909      _node_num = countNodes(_graph);
910      _arc_num = countArcs(_graph);
911      int all_node_num = _node_num + 1;
912      int max_arc_num = _arc_num + 2 * _node_num;
913
914      _source.resize(max_arc_num);
915      _target.resize(max_arc_num);
916
917      _lower.resize(_arc_num);
918      _upper.resize(_arc_num);
919      _cap.resize(max_arc_num);
920      _cost.resize(max_arc_num);
921      _supply.resize(all_node_num);
922      _flow.resize(max_arc_num);
923      _pi.resize(all_node_num);
924
925      _parent.resize(all_node_num);
926      _pred.resize(all_node_num);
[990]927      _pred_dir.resize(all_node_num);
[898]928      _thread.resize(all_node_num);
929      _rev_thread.resize(all_node_num);
930      _succ_num.resize(all_node_num);
931      _last_succ.resize(all_node_num);
932      _state.resize(max_arc_num);
933
934      // Copy the graph
935      int i = 0;
936      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
937        _node_id[n] = i;
938      }
939      if (_arc_mixing) {
940        // Store the arcs in a mixed order
[991]941        const int skip = std::max(_arc_num / _node_num, 3);
[898]942        int i = 0, j = 0;
943        for (ArcIt a(_graph); a != INVALID; ++a) {
944          _arc_id[a] = i;
945          _source[i] = _node_id[_graph.source(a)];
946          _target[i] = _node_id[_graph.target(a)];
[991]947          if ((i += skip) >= _arc_num) i = ++j;
[898]948        }
949      } else {
950        // Store the arcs in the original order
951        int i = 0;
952        for (ArcIt a(_graph); a != INVALID; ++a, ++i) {
953          _arc_id[a] = i;
954          _source[i] = _node_id[_graph.source(a)];
955          _target[i] = _node_id[_graph.target(a)];
956        }
957      }
[956]958
[898]959      // Reset parameters
960      resetParams();
961      return *this;
962    }
[956]963
[648]964    /// @}
965
966    /// \name Query Functions
967    /// The results of the algorithm can be obtained using these
968    /// functions.\n
[652]969    /// The \ref run() function must be called before using them.
970
[648]971    /// @{
972
[652]973    /// \brief Return the total cost of the found flow.
974    ///
975    /// This function returns the total cost of the found flow.
[1254]976    /// Its complexity is O(m).
[652]977    ///
978    /// \note The return type of the function can be specified as a
979    /// template parameter. For example,
980    /// \code
981    ///   ns.totalCost<double>();
982    /// \endcode
[654]983    /// It is useful if the total cost cannot be stored in the \c Cost
[652]984    /// type of the algorithm, which is the default return type of the
985    /// function.
986    ///
987    /// \pre \ref run() must be called before using this function.
[689]988    template <typename Number>
989    Number totalCost() const {
990      Number c = 0;
991      for (ArcIt a(_graph); a != INVALID; ++a) {
992        int i = _arc_id[a];
993        c += Number(_flow[i]) * Number(_cost[i]);
[652]994      }
995      return c;
996    }
997
998#ifndef DOXYGEN
[654]999    Cost totalCost() const {
1000      return totalCost<Cost>();
[652]1001    }
1002#endif
1003
1004    /// \brief Return the flow on the given arc.
1005    ///
1006    /// This function returns the flow on the given arc.
1007    ///
1008    /// \pre \ref run() must be called before using this function.
[688]1009    Value flow(const Arc& a) const {
[689]1010      return _flow[_arc_id[a]];
[652]1011    }
1012
[1165]1013    /// \brief Copy the flow values (the primal solution) into the
1014    /// given map.
[648]1015    ///
[689]1016    /// This function copies the flow value on each arc into the given
1017    /// map. The \c Value type of the algorithm must be convertible to
1018    /// the \c Value type of the map.
[648]1019    ///
1020    /// \pre \ref run() must be called before using this function.
[689]1021    template <typename FlowMap>
1022    void flowMap(FlowMap &map) const {
1023      for (ArcIt a(_graph); a != INVALID; ++a) {
1024        map.set(a, _flow[_arc_id[a]]);
1025      }
[648]1026    }
1027
[652]1028    /// \brief Return the potential (dual value) of the given node.
1029    ///
1030    /// This function returns the potential (dual value) of the
1031    /// given node.
1032    ///
1033    /// \pre \ref run() must be called before using this function.
[654]1034    Cost potential(const Node& n) const {
[689]1035      return _pi[_node_id[n]];
[652]1036    }
1037
[1165]1038    /// \brief Copy the potential values (the dual solution) into the
1039    /// given map.
[648]1040    ///
[689]1041    /// This function copies the potential (dual value) of each node
1042    /// into the given map.
1043    /// The \c Cost type of the algorithm must be convertible to the
1044    /// \c Value type of the map.
[648]1045    ///
1046    /// \pre \ref run() must be called before using this function.
[689]1047    template <typename PotentialMap>
1048    void potentialMap(PotentialMap &map) const {
1049      for (NodeIt n(_graph); n != INVALID; ++n) {
1050        map.set(n, _pi[_node_id[n]]);
1051      }
[648]1052    }
1053
1054    /// @}
1055
1056  private:
1057
1058    // Initialize internal data structures
1059    bool init() {
[652]1060      if (_node_num == 0) return false;
[648]1061
[689]1062      // Check the sum of supply values
1063      _sum_supply = 0;
1064      for (int i = 0; i != _node_num; ++i) {
1065        _sum_supply += _supply[i];
1066      }
[690]1067      if ( !((_stype == GEQ && _sum_supply <= 0) ||
1068             (_stype == LEQ && _sum_supply >= 0)) ) return false;
[648]1069
[1240]1070      // Check lower and upper bounds
1071      LEMON_DEBUG(checkBoundMaps(),
1072          "Upper bounds must be greater or equal to the lower bounds");
1073
[689]1074      // Remove non-zero lower bounds
1075      if (_have_lower) {
1076        for (int i = 0; i != _arc_num; ++i) {
1077          Value c = _lower[i];
1078          if (c >= 0) {
[877]1079            _cap[i] = _upper[i] < MAX ? _upper[i] - c : INF;
[689]1080          } else {
[877]1081            _cap[i] = _upper[i] < MAX + c ? _upper[i] - c : INF;
[689]1082          }
1083          _supply[_source[i]] -= c;
1084          _supply[_target[i]] += c;
1085        }
1086      } else {
1087        for (int i = 0; i != _arc_num; ++i) {
1088          _cap[i] = _upper[i];
1089        }
[652]1090      }
[648]1091
[656]1092      // Initialize artifical cost
[687]1093      Cost ART_COST;
[656]1094      if (std::numeric_limits<Cost>::is_exact) {
[710]1095        ART_COST = std::numeric_limits<Cost>::max() / 2 + 1;
[656]1096      } else {
[976]1097        ART_COST = 0;
[656]1098        for (int i = 0; i != _arc_num; ++i) {
[687]1099          if (_cost[i] > ART_COST) ART_COST = _cost[i];
[656]1100        }
[687]1101        ART_COST = (ART_COST + 1) * _node_num;
[656]1102      }
1103
[689]1104      // Initialize arc maps
1105      for (int i = 0; i != _arc_num; ++i) {
1106        _flow[i] = 0;
1107        _state[i] = STATE_LOWER;
1108      }
[956]1109
[648]1110      // Set data for the artificial root node
1111      _root = _node_num;
1112      _parent[_root] = -1;
1113      _pred[_root] = -1;
1114      _thread[_root] = 0;
[651]1115      _rev_thread[0] = _root;
[689]1116      _succ_num[_root] = _node_num + 1;
[651]1117      _last_succ[_root] = _root - 1;
[687]1118      _supply[_root] = -_sum_supply;
[710]1119      _pi[_root] = 0;
[648]1120
1121      // Add artificial arcs and initialize the spanning tree data structure
[710]1122      if (_sum_supply == 0) {
1123        // EQ supply constraints
1124        _search_arc_num = _arc_num;
1125        _all_arc_num = _arc_num + _node_num;
1126        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1127          _parent[u] = _root;
1128          _pred[u] = e;
1129          _thread[u] = u + 1;
1130          _rev_thread[u + 1] = u;
1131          _succ_num[u] = 1;
1132          _last_succ[u] = u;
1133          _cap[e] = INF;
1134          _state[e] = STATE_TREE;
1135          if (_supply[u] >= 0) {
[990]1136            _pred_dir[u] = DIR_UP;
[710]1137            _pi[u] = 0;
1138            _source[e] = u;
1139            _target[e] = _root;
1140            _flow[e] = _supply[u];
1141            _cost[e] = 0;
1142          } else {
[990]1143            _pred_dir[u] = DIR_DOWN;
[710]1144            _pi[u] = ART_COST;
1145            _source[e] = _root;
1146            _target[e] = u;
1147            _flow[e] = -_supply[u];
1148            _cost[e] = ART_COST;
1149          }
[648]1150        }
1151      }
[710]1152      else if (_sum_supply > 0) {
1153        // LEQ supply constraints
1154        _search_arc_num = _arc_num + _node_num;
1155        int f = _arc_num + _node_num;
1156        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1157          _parent[u] = _root;
1158          _thread[u] = u + 1;
1159          _rev_thread[u + 1] = u;
1160          _succ_num[u] = 1;
1161          _last_succ[u] = u;
1162          if (_supply[u] >= 0) {
[990]1163            _pred_dir[u] = DIR_UP;
[710]1164            _pi[u] = 0;
1165            _pred[u] = e;
1166            _source[e] = u;
1167            _target[e] = _root;
1168            _cap[e] = INF;
1169            _flow[e] = _supply[u];
1170            _cost[e] = 0;
1171            _state[e] = STATE_TREE;
1172          } else {
[990]1173            _pred_dir[u] = DIR_DOWN;
[710]1174            _pi[u] = ART_COST;
1175            _pred[u] = f;
1176            _source[f] = _root;
1177            _target[f] = u;
1178            _cap[f] = INF;
1179            _flow[f] = -_supply[u];
1180            _cost[f] = ART_COST;
1181            _state[f] = STATE_TREE;
1182            _source[e] = u;
1183            _target[e] = _root;
1184            _cap[e] = INF;
1185            _flow[e] = 0;
1186            _cost[e] = 0;
1187            _state[e] = STATE_LOWER;
1188            ++f;
1189          }
1190        }
1191        _all_arc_num = f;
1192      }
1193      else {
1194        // GEQ supply constraints
1195        _search_arc_num = _arc_num + _node_num;
1196        int f = _arc_num + _node_num;
1197        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1198          _parent[u] = _root;
1199          _thread[u] = u + 1;
1200          _rev_thread[u + 1] = u;
1201          _succ_num[u] = 1;
1202          _last_succ[u] = u;
1203          if (_supply[u] <= 0) {
[990]1204            _pred_dir[u] = DIR_DOWN;
[710]1205            _pi[u] = 0;
1206            _pred[u] = e;
1207            _source[e] = _root;
1208            _target[e] = u;
1209            _cap[e] = INF;
1210            _flow[e] = -_supply[u];
1211            _cost[e] = 0;
1212            _state[e] = STATE_TREE;
1213          } else {
[990]1214            _pred_dir[u] = DIR_UP;
[710]1215            _pi[u] = -ART_COST;
1216            _pred[u] = f;
1217            _source[f] = u;
1218            _target[f] = _root;
1219            _cap[f] = INF;
1220            _flow[f] = _supply[u];
1221            _state[f] = STATE_TREE;
1222            _cost[f] = ART_COST;
1223            _source[e] = _root;
1224            _target[e] = u;
1225            _cap[e] = INF;
1226            _flow[e] = 0;
1227            _cost[e] = 0;
1228            _state[e] = STATE_LOWER;
1229            ++f;
1230          }
1231        }
1232        _all_arc_num = f;
1233      }
[648]1234
1235      return true;
1236    }
[1270]1237
[1240]1238    // Check if the upper bound is greater or equal to the lower bound
1239    // on each arc.
1240    bool checkBoundMaps() {
1241      for (int j = 0; j != _arc_num; ++j) {
1242        if (_upper[j] < _lower[j]) return false;
1243      }
1244      return true;
1245    }
[648]1246
1247    // Find the join node
1248    void findJoinNode() {
[650]1249      int u = _source[in_arc];
1250      int v = _target[in_arc];
[648]1251      while (u != v) {
[651]1252        if (_succ_num[u] < _succ_num[v]) {
1253          u = _parent[u];
1254        } else {
1255          v = _parent[v];
1256        }
[648]1257      }
1258      join = u;
1259    }
1260
1261    // Find the leaving arc of the cycle and returns true if the
1262    // leaving arc is not the same as the entering arc
1263    bool findLeavingArc() {
1264      // Initialize first and second nodes according to the direction
1265      // of the cycle
[990]1266      int first, second;
[650]1267      if (_state[in_arc] == STATE_LOWER) {
1268        first  = _source[in_arc];
1269        second = _target[in_arc];
[648]1270      } else {
[650]1271        first  = _target[in_arc];
1272        second = _source[in_arc];
[648]1273      }
[650]1274      delta = _cap[in_arc];
[648]1275      int result = 0;
[990]1276      Value c, d;
[648]1277      int e;
1278
[990]1279      // Search the cycle form the first node to the join node
[648]1280      for (int u = first; u != join; u = _parent[u]) {
1281        e = _pred[u];
[990]1282        d = _flow[e];
1283        if (_pred_dir[u] == DIR_DOWN) {
1284          c = _cap[e];
1285          d = c >= MAX ? INF : c - d;
1286        }
[648]1287        if (d < delta) {
1288          delta = d;
1289          u_out = u;
1290          result = 1;
1291        }
1292      }
[990]1293
1294      // Search the cycle form the second node to the join node
[648]1295      for (int u = second; u != join; u = _parent[u]) {
1296        e = _pred[u];
[990]1297        d = _flow[e];
1298        if (_pred_dir[u] == DIR_UP) {
1299          c = _cap[e];
1300          d = c >= MAX ? INF : c - d;
1301        }
[648]1302        if (d <= delta) {
1303          delta = d;
1304          u_out = u;
1305          result = 2;
1306        }
1307      }
1308
1309      if (result == 1) {
1310        u_in = first;
1311        v_in = second;
1312      } else {
1313        u_in = second;
1314        v_in = first;
1315      }
1316      return result != 0;
1317    }
1318
1319    // Change _flow and _state vectors
1320    void changeFlow(bool change) {
1321      // Augment along the cycle
1322      if (delta > 0) {
[688]1323        Value val = _state[in_arc] * delta;
[650]1324        _flow[in_arc] += val;
1325        for (int u = _source[in_arc]; u != join; u = _parent[u]) {
[990]1326          _flow[_pred[u]] -= _pred_dir[u] * val;
[648]1327        }
[650]1328        for (int u = _target[in_arc]; u != join; u = _parent[u]) {
[990]1329          _flow[_pred[u]] += _pred_dir[u] * val;
[648]1330        }
1331      }
1332      // Update the state of the entering and leaving arcs
1333      if (change) {
[650]1334        _state[in_arc] = STATE_TREE;
[648]1335        _state[_pred[u_out]] =
1336          (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
1337      } else {
[650]1338        _state[in_arc] = -_state[in_arc];
[648]1339      }
1340    }
1341
[651]1342    // Update the tree structure
1343    void updateTreeStructure() {
1344      int old_rev_thread = _rev_thread[u_out];
1345      int old_succ_num = _succ_num[u_out];
1346      int old_last_succ = _last_succ[u_out];
[648]1347      v_out = _parent[u_out];
1348
[990]1349      // Check if u_in and u_out coincide
1350      if (u_in == u_out) {
1351        // Update _parent, _pred, _pred_dir
1352        _parent[u_in] = v_in;
1353        _pred[u_in] = in_arc;
1354        _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
[651]1355
[990]1356        // Update _thread and _rev_thread
1357        if (_thread[v_in] != u_out) {
1358          int after = _thread[old_last_succ];
1359          _thread[old_rev_thread] = after;
1360          _rev_thread[after] = old_rev_thread;
1361          after = _thread[v_in];
1362          _thread[v_in] = u_out;
1363          _rev_thread[u_out] = v_in;
1364          _thread[old_last_succ] = after;
1365          _rev_thread[after] = old_last_succ;
1366        }
[651]1367      } else {
[990]1368        // Handle the case when old_rev_thread equals to v_in
1369        // (it also means that join and v_out coincide)
1370        int thread_continue = old_rev_thread == v_in ?
1371          _thread[old_last_succ] : _thread[v_in];
[648]1372
[990]1373        // Update _thread and _parent along the stem nodes (i.e. the nodes
1374        // between u_in and u_out, whose parent have to be changed)
1375        int stem = u_in;              // the current stem node
1376        int par_stem = v_in;          // the new parent of stem
1377        int next_stem;                // the next stem node
1378        int last = _last_succ[u_in];  // the last successor of stem
1379        int before, after = _thread[last];
1380        _thread[v_in] = u_in;
1381        _dirty_revs.clear();
1382        _dirty_revs.push_back(v_in);
1383        while (stem != u_out) {
1384          // Insert the next stem node into the thread list
1385          next_stem = _parent[stem];
1386          _thread[last] = next_stem;
1387          _dirty_revs.push_back(last);
[648]1388
[990]1389          // Remove the subtree of stem from the thread list
1390          before = _rev_thread[stem];
1391          _thread[before] = after;
1392          _rev_thread[after] = before;
[648]1393
[990]1394          // Change the parent node and shift stem nodes
1395          _parent[stem] = par_stem;
1396          par_stem = stem;
1397          stem = next_stem;
[648]1398
[990]1399          // Update last and after
1400          last = _last_succ[stem] == _last_succ[par_stem] ?
1401            _rev_thread[par_stem] : _last_succ[stem];
1402          after = _thread[last];
1403        }
1404        _parent[u_out] = par_stem;
1405        _thread[last] = thread_continue;
1406        _rev_thread[thread_continue] = last;
1407        _last_succ[u_out] = last;
[648]1408
[990]1409        // Remove the subtree of u_out from the thread list except for
1410        // the case when old_rev_thread equals to v_in
1411        if (old_rev_thread != v_in) {
1412          _thread[old_rev_thread] = after;
1413          _rev_thread[after] = old_rev_thread;
1414        }
[651]1415
[990]1416        // Update _rev_thread using the new _thread values
1417        for (int i = 0; i != int(_dirty_revs.size()); ++i) {
1418          int u = _dirty_revs[i];
1419          _rev_thread[_thread[u]] = u;
1420        }
[651]1421
[990]1422        // Update _pred, _pred_dir, _last_succ and _succ_num for the
1423        // stem nodes from u_out to u_in
1424        int tmp_sc = 0, tmp_ls = _last_succ[u_out];
1425        for (int u = u_out, p = _parent[u]; u != u_in; u = p, p = _parent[u]) {
1426          _pred[u] = _pred[p];
1427          _pred_dir[u] = -_pred_dir[p];
1428          tmp_sc += _succ_num[u] - _succ_num[p];
1429          _succ_num[u] = tmp_sc;
1430          _last_succ[p] = tmp_ls;
1431        }
1432        _pred[u_in] = in_arc;
1433        _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
1434        _succ_num[u_in] = old_succ_num;
[651]1435      }
1436
1437      // Update _last_succ from v_in towards the root
[990]1438      int up_limit_out = _last_succ[join] == v_in ? join : -1;
1439      int last_succ_out = _last_succ[u_out];
1440      for (int u = v_in; u != -1 && _last_succ[u] == v_in; u = _parent[u]) {
1441        _last_succ[u] = last_succ_out;
[651]1442      }
[990]1443
[651]1444      // Update _last_succ from v_out towards the root
1445      if (join != old_rev_thread && v_in != old_rev_thread) {
[990]1446        for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
[651]1447             u = _parent[u]) {
1448          _last_succ[u] = old_rev_thread;
1449        }
[990]1450      }
1451      else if (last_succ_out != old_last_succ) {
1452        for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
[651]1453             u = _parent[u]) {
[990]1454          _last_succ[u] = last_succ_out;
[651]1455        }
1456      }
1457
1458      // Update _succ_num from v_in to join
[990]1459      for (int u = v_in; u != join; u = _parent[u]) {
[651]1460        _succ_num[u] += old_succ_num;
1461      }
1462      // Update _succ_num from v_out to join
[990]1463      for (int u = v_out; u != join; u = _parent[u]) {
[651]1464        _succ_num[u] -= old_succ_num;
[648]1465      }
1466    }
1467
[990]1468    // Update potentials in the subtree that has been moved
[651]1469    void updatePotential() {
[990]1470      Cost sigma = _pi[v_in] - _pi[u_in] -
1471                   _pred_dir[u_in] * _cost[in_arc];
[655]1472      int end = _thread[_last_succ[u_in]];
1473      for (int u = u_in; u != end; u = _thread[u]) {
1474        _pi[u] += sigma;
[648]1475      }
1476    }
1477
[910]1478    // Heuristic initial pivots
1479    bool initialPivots() {
1480      Value curr, total = 0;
1481      std::vector<Node> supply_nodes, demand_nodes;
1482      for (NodeIt u(_graph); u != INVALID; ++u) {
1483        curr = _supply[_node_id[u]];
1484        if (curr > 0) {
1485          total += curr;
1486          supply_nodes.push_back(u);
1487        }
1488        else if (curr < 0) {
1489          demand_nodes.push_back(u);
1490        }
1491      }
1492      if (_sum_supply > 0) total -= _sum_supply;
1493      if (total <= 0) return true;
1494
1495      IntVector arc_vector;
1496      if (_sum_supply >= 0) {
1497        if (supply_nodes.size() == 1 && demand_nodes.size() == 1) {
1498          // Perform a reverse graph search from the sink to the source
1499          typename GR::template NodeMap<bool> reached(_graph, false);
1500          Node s = supply_nodes[0], t = demand_nodes[0];
1501          std::vector<Node> stack;
1502          reached[t] = true;
1503          stack.push_back(t);
1504          while (!stack.empty()) {
1505            Node u, v = stack.back();
1506            stack.pop_back();
1507            if (v == s) break;
1508            for (InArcIt a(_graph, v); a != INVALID; ++a) {
1509              if (reached[u = _graph.source(a)]) continue;
1510              int j = _arc_id[a];
1511              if (_cap[j] >= total) {
1512                arc_vector.push_back(j);
1513                reached[u] = true;
1514                stack.push_back(u);
1515              }
1516            }
1517          }
1518        } else {
[1217]1519          // Find the min. cost incoming arc for each demand node
[910]1520          for (int i = 0; i != int(demand_nodes.size()); ++i) {
1521            Node v = demand_nodes[i];
1522            Cost c, min_cost = std::numeric_limits<Cost>::max();
1523            Arc min_arc = INVALID;
1524            for (InArcIt a(_graph, v); a != INVALID; ++a) {
1525              c = _cost[_arc_id[a]];
1526              if (c < min_cost) {
1527                min_cost = c;
1528                min_arc = a;
1529              }
1530            }
1531            if (min_arc != INVALID) {
1532              arc_vector.push_back(_arc_id[min_arc]);
1533            }
1534          }
1535        }
1536      } else {
1537        // Find the min. cost outgoing arc for each supply node
1538        for (int i = 0; i != int(supply_nodes.size()); ++i) {
1539          Node u = supply_nodes[i];
1540          Cost c, min_cost = std::numeric_limits<Cost>::max();
1541          Arc min_arc = INVALID;
1542          for (OutArcIt a(_graph, u); a != INVALID; ++a) {
1543            c = _cost[_arc_id[a]];
1544            if (c < min_cost) {
1545              min_cost = c;
1546              min_arc = a;
1547            }
1548          }
1549          if (min_arc != INVALID) {
1550            arc_vector.push_back(_arc_id[min_arc]);
1551          }
1552        }
1553      }
1554
1555      // Perform heuristic initial pivots
1556      for (int i = 0; i != int(arc_vector.size()); ++i) {
1557        in_arc = arc_vector[i];
1558        if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] -
1559            _pi[_target[in_arc]]) >= 0) continue;
1560        findJoinNode();
1561        bool change = findLeavingArc();
1562        if (delta >= MAX) return false;
1563        changeFlow(change);
1564        if (change) {
1565          updateTreeStructure();
1566          updatePotential();
1567        }
1568      }
1569      return true;
1570    }
1571
[648]1572    // Execute the algorithm
[687]1573    ProblemType start(PivotRule pivot_rule) {
[648]1574      // Select the pivot rule implementation
1575      switch (pivot_rule) {
[652]1576        case FIRST_ELIGIBLE:
[648]1577          return start<FirstEligiblePivotRule>();
[652]1578        case BEST_ELIGIBLE:
[648]1579          return start<BestEligiblePivotRule>();
[652]1580        case BLOCK_SEARCH:
[648]1581          return start<BlockSearchPivotRule>();
[652]1582        case CANDIDATE_LIST:
[648]1583          return start<CandidateListPivotRule>();
[652]1584        case ALTERING_LIST:
[648]1585          return start<AlteringListPivotRule>();
1586      }
[687]1587      return INFEASIBLE; // avoid warning
[648]1588    }
1589
[652]1590    template <typename PivotRuleImpl>
[687]1591    ProblemType start() {
[652]1592      PivotRuleImpl pivot(*this);
[648]1593
[910]1594      // Perform heuristic initial pivots
1595      if (!initialPivots()) return UNBOUNDED;
1596
[652]1597      // Execute the Network Simplex algorithm
[648]1598      while (pivot.findEnteringArc()) {
1599        findJoinNode();
1600        bool change = findLeavingArc();
[877]1601        if (delta >= MAX) return UNBOUNDED;
[648]1602        changeFlow(change);
1603        if (change) {
[651]1604          updateTreeStructure();
1605          updatePotential();
[648]1606        }
1607      }
[956]1608
[687]1609      // Check feasibility
[710]1610      for (int e = _search_arc_num; e != _all_arc_num; ++e) {
1611        if (_flow[e] != 0) return INFEASIBLE;
[687]1612      }
[648]1613
[689]1614      // Transform the solution and the supply map to the original form
1615      if (_have_lower) {
[648]1616        for (int i = 0; i != _arc_num; ++i) {
[689]1617          Value c = _lower[i];
1618          if (c != 0) {
1619            _flow[i] += c;
1620            _supply[_source[i]] += c;
1621            _supply[_target[i]] -= c;
1622          }
[648]1623        }
1624      }
[956]1625
[710]1626      // Shift potentials to meet the requirements of the GEQ/LEQ type
1627      // optimality conditions
1628      if (_sum_supply == 0) {
1629        if (_stype == GEQ) {
[976]1630          Cost max_pot = -std::numeric_limits<Cost>::max();
[710]1631          for (int i = 0; i != _node_num; ++i) {
1632            if (_pi[i] > max_pot) max_pot = _pi[i];
1633          }
1634          if (max_pot > 0) {
1635            for (int i = 0; i != _node_num; ++i)
1636              _pi[i] -= max_pot;
1637          }
1638        } else {
1639          Cost min_pot = std::numeric_limits<Cost>::max();
1640          for (int i = 0; i != _node_num; ++i) {
1641            if (_pi[i] < min_pot) min_pot = _pi[i];
1642          }
1643          if (min_pot < 0) {
1644            for (int i = 0; i != _node_num; ++i)
1645              _pi[i] -= min_pot;
1646          }
1647        }
1648      }
[648]1649
[687]1650      return OPTIMAL;
[648]1651    }
1652
1653  }; //class NetworkSimplex
1654
1655  ///@}
1656
1657} //namespace lemon
1658
1659#endif //LEMON_NETWORK_SIMPLEX_H
Note: See TracBrowser for help on using the repository browser.