COIN-OR::LEMON - Graph Library

source: lemon/lemon/network_simplex.h @ 1136:fcb6ad1e67d0

Last change on this file since 1136:fcb6ad1e67d0 was 1136:fcb6ad1e67d0, checked in by Peter Kovacs <kpeter@…>, 8 years ago

Improve the Altering List pivot rule for NetworkSimplex? (#435)
Much less candidate arcs are preserved from an iteration to the
next one and partial_sort() is used instead of heap operations.

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