COIN-OR::LEMON - Graph Library

source: lemon/lemon/network_simplex.h @ 656:e6927fe719e6

Last change on this file since 656:e6927fe719e6 was 656:e6927fe719e6, checked in by Peter Kovacs <kpeter@…>, 15 years ago

Support >= and <= constraints in NetworkSimplex? (#219, #234)

By default the same inequality constraints are supported as by
Circulation (the GEQ form), but the LEQ form can also be selected
using the problemType() function.

The documentation of the min. cost flow module is reworked and
extended with important notes and explanations about the different
variants of the problem and about the dual solution and optimality
conditions.

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