COIN-OR::LEMON - Graph Library

source: lemon/lemon/network_simplex.h @ 1165:16f55008c863

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

Doc improvements for min cost flow algorithms (#437)

File size: 51.0 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 solving this problem.
52  /// (For more information, see \ref min_cost_flow_algs "the module page".)
53  /// Furthermore, this class supports both directions of the supply/demand
54  /// inequality constraints. For more information, see \ref SupplyType.
55  ///
56  /// Most of the parameters of the problem (except for the digraph)
57  /// can be given using separate functions, and the algorithm can be
58  /// executed using the \ref run() function. If some parameters are not
59  /// specified, then default values will be used.
60  ///
61  /// \tparam GR The digraph type the algorithm runs on.
62  /// \tparam V The number type used for flow amounts, capacity bounds
63  /// and supply values in the algorithm. By default, it is \c int.
64  /// \tparam C The number type used for costs and potentials in the
65  /// algorithm. By default, it is the same as \c V.
66  ///
67  /// \warning Both \c V and \c C must be signed number types.
68  /// \warning All input data (capacities, supply values, and costs) must
69  /// be integer.
70  ///
71  /// \note %NetworkSimplex provides five different pivot rule
72  /// implementations, from which the most efficient one is used
73  /// by default. For more information, see \ref PivotRule.
74  template <typename GR, typename V = int, typename C = V>
75  class NetworkSimplex
76  {
77  public:
78
79    /// The type of the flow amounts, capacity bounds and supply values
80    typedef V Value;
81    /// The type of the arc costs
82    typedef C Cost;
83
84  public:
85
86    /// \brief Problem type constants for the \c run() function.
87    ///
88    /// Enum type containing the problem type constants that can be
89    /// returned by the \ref run() function of the algorithm.
90    enum ProblemType {
91      /// The problem has no feasible solution (flow).
92      INFEASIBLE,
93      /// The problem has optimal solution (i.e. it is feasible and
94      /// bounded), and the algorithm has found optimal flow and node
95      /// potentials (primal and dual solutions).
96      OPTIMAL,
97      /// The objective function of the problem is unbounded, i.e.
98      /// there is a directed cycle having negative total cost and
99      /// infinite upper bound.
100      UNBOUNDED
101    };
102
103    /// \brief Constants for selecting the type of the supply constraints.
104    ///
105    /// Enum type containing constants for selecting the supply type,
106    /// i.e. the direction of the inequalities in the supply/demand
107    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
108    ///
109    /// The default supply type is \c GEQ, the \c LEQ type can be
110    /// selected using \ref supplyType().
111    /// The equality form is a special case of both supply types.
112    enum SupplyType {
113      /// This option means that there are <em>"greater or equal"</em>
114      /// supply/demand constraints in the definition of the problem.
115      GEQ,
116      /// This option means that there are <em>"less or equal"</em>
117      /// supply/demand constraints in the definition of the problem.
118      LEQ
119    };
120
121    /// \brief Constants for selecting the pivot rule.
122    ///
123    /// Enum type containing constants for selecting the pivot rule for
124    /// the \ref run() function.
125    ///
126    /// \ref NetworkSimplex provides five different pivot rule
127    /// implementations that significantly affect the running time
128    /// of the algorithm.
129    /// By default, \ref BLOCK_SEARCH "Block Search" is used, which
130    /// turend out to be the most efficient and the most robust on various
131    /// test inputs.
132    /// However, another pivot rule can be selected using the \ref run()
133    /// function with the proper parameter.
134    enum PivotRule {
135
136      /// The \e First \e Eligible pivot rule.
137      /// The next eligible arc is selected in a wraparound fashion
138      /// in every iteration.
139      FIRST_ELIGIBLE,
140
141      /// The \e Best \e Eligible pivot rule.
142      /// The best eligible arc is selected in every iteration.
143      BEST_ELIGIBLE,
144
145      /// The \e Block \e Search pivot rule.
146      /// A specified number of arcs are examined in every iteration
147      /// in a wraparound fashion and the best eligible arc is selected
148      /// from this block.
149      BLOCK_SEARCH,
150
151      /// The \e Candidate \e List pivot rule.
152      /// In a major iteration a candidate list is built from eligible arcs
153      /// in a wraparound fashion and in the following minor iterations
154      /// the best eligible arc is selected from this list.
155      CANDIDATE_LIST,
156
157      /// The \e Altering \e Candidate \e List pivot rule.
158      /// It is a modified version of the Candidate List method.
159      /// It keeps only the several best eligible arcs from the former
160      /// candidate list and extends this list in every iteration.
161      ALTERING_LIST
162    };
163
164  private:
165
166    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
167
168    typedef std::vector<int> IntVector;
169    typedef std::vector<Value> ValueVector;
170    typedef std::vector<Cost> CostVector;
171    typedef std::vector<signed char> CharVector;
172    // Note: vector<signed char> is used instead of vector<ArcState> and
173    // vector<ArcDirection> for efficiency reasons
174
175    // State constants for arcs
176    enum ArcState {
177      STATE_UPPER = -1,
178      STATE_TREE  =  0,
179      STATE_LOWER =  1
180    };
181
182    // Direction constants for tree arcs
183    enum ArcDirection {
184      DIR_DOWN = -1,
185      DIR_UP   =  1
186    };
187
188  private:
189
190    // Data related to the underlying digraph
191    const GR &_graph;
192    int _node_num;
193    int _arc_num;
194    int _all_arc_num;
195    int _search_arc_num;
196
197    // Parameters of the problem
198    bool _have_lower;
199    SupplyType _stype;
200    Value _sum_supply;
201
202    // Data structures for storing the digraph
203    IntNodeMap _node_id;
204    IntArcMap _arc_id;
205    IntVector _source;
206    IntVector _target;
207    bool _arc_mixing;
208
209    // Node and arc data
210    ValueVector _lower;
211    ValueVector _upper;
212    ValueVector _cap;
213    CostVector _cost;
214    ValueVector _supply;
215    ValueVector _flow;
216    CostVector _pi;
217
218    // Data for storing the spanning tree structure
219    IntVector _parent;
220    IntVector _pred;
221    IntVector _thread;
222    IntVector _rev_thread;
223    IntVector _succ_num;
224    IntVector _last_succ;
225    CharVector _pred_dir;
226    CharVector _state;
227    IntVector _dirty_revs;
228    int _root;
229
230    // Temporary data used in the current pivot iteration
231    int in_arc, join, u_in, v_in, u_out, v_out;
232    Value delta;
233
234    const Value MAX;
235
236  public:
237
238    /// \brief Constant for infinite upper bounds (capacities).
239    ///
240    /// Constant for infinite upper bounds (capacities).
241    /// It is \c std::numeric_limits<Value>::infinity() if available,
242    /// \c std::numeric_limits<Value>::max() otherwise.
243    const Value INF;
244
245  private:
246
247    // Implementation of the First Eligible pivot rule
248    class FirstEligiblePivotRule
249    {
250    private:
251
252      // References to the NetworkSimplex class
253      const IntVector  &_source;
254      const IntVector  &_target;
255      const CostVector &_cost;
256      const CharVector &_state;
257      const CostVector &_pi;
258      int &_in_arc;
259      int _search_arc_num;
260
261      // Pivot rule data
262      int _next_arc;
263
264    public:
265
266      // Constructor
267      FirstEligiblePivotRule(NetworkSimplex &ns) :
268        _source(ns._source), _target(ns._target),
269        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
270        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
271        _next_arc(0)
272      {}
273
274      // Find next entering arc
275      bool findEnteringArc() {
276        Cost c;
277        for (int e = _next_arc; e != _search_arc_num; ++e) {
278          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
279          if (c < 0) {
280            _in_arc = e;
281            _next_arc = e + 1;
282            return true;
283          }
284        }
285        for (int e = 0; e != _next_arc; ++e) {
286          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
287          if (c < 0) {
288            _in_arc = e;
289            _next_arc = e + 1;
290            return true;
291          }
292        }
293        return false;
294      }
295
296    }; //class FirstEligiblePivotRule
297
298
299    // Implementation of the Best Eligible pivot rule
300    class BestEligiblePivotRule
301    {
302    private:
303
304      // References to the NetworkSimplex class
305      const IntVector  &_source;
306      const IntVector  &_target;
307      const CostVector &_cost;
308      const CharVector &_state;
309      const CostVector &_pi;
310      int &_in_arc;
311      int _search_arc_num;
312
313    public:
314
315      // Constructor
316      BestEligiblePivotRule(NetworkSimplex &ns) :
317        _source(ns._source), _target(ns._target),
318        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
319        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num)
320      {}
321
322      // Find next entering arc
323      bool findEnteringArc() {
324        Cost c, min = 0;
325        for (int e = 0; e != _search_arc_num; ++e) {
326          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
327          if (c < min) {
328            min = c;
329            _in_arc = e;
330          }
331        }
332        return min < 0;
333      }
334
335    }; //class BestEligiblePivotRule
336
337
338    // Implementation of the Block Search pivot rule
339    class BlockSearchPivotRule
340    {
341    private:
342
343      // References to the NetworkSimplex class
344      const IntVector  &_source;
345      const IntVector  &_target;
346      const CostVector &_cost;
347      const CharVector &_state;
348      const CostVector &_pi;
349      int &_in_arc;
350      int _search_arc_num;
351
352      // Pivot rule data
353      int _block_size;
354      int _next_arc;
355
356    public:
357
358      // Constructor
359      BlockSearchPivotRule(NetworkSimplex &ns) :
360        _source(ns._source), _target(ns._target),
361        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
362        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
363        _next_arc(0)
364      {
365        // The main parameters of the pivot rule
366        const double BLOCK_SIZE_FACTOR = 1.0;
367        const int MIN_BLOCK_SIZE = 10;
368
369        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
370                                    std::sqrt(double(_search_arc_num))),
371                                MIN_BLOCK_SIZE );
372      }
373
374      // Find next entering arc
375      bool findEnteringArc() {
376        Cost c, min = 0;
377        int cnt = _block_size;
378        int e;
379        for (e = _next_arc; e != _search_arc_num; ++e) {
380          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
381          if (c < min) {
382            min = c;
383            _in_arc = e;
384          }
385          if (--cnt == 0) {
386            if (min < 0) goto search_end;
387            cnt = _block_size;
388          }
389        }
390        for (e = 0; e != _next_arc; ++e) {
391          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
392          if (c < min) {
393            min = c;
394            _in_arc = e;
395          }
396          if (--cnt == 0) {
397            if (min < 0) goto search_end;
398            cnt = _block_size;
399          }
400        }
401        if (min >= 0) return false;
402
403      search_end:
404        _next_arc = e;
405        return true;
406      }
407
408    }; //class BlockSearchPivotRule
409
410
411    // Implementation of the Candidate List pivot rule
412    class CandidateListPivotRule
413    {
414    private:
415
416      // References to the NetworkSimplex class
417      const IntVector  &_source;
418      const IntVector  &_target;
419      const CostVector &_cost;
420      const CharVector &_state;
421      const CostVector &_pi;
422      int &_in_arc;
423      int _search_arc_num;
424
425      // Pivot rule data
426      IntVector _candidates;
427      int _list_length, _minor_limit;
428      int _curr_length, _minor_count;
429      int _next_arc;
430
431    public:
432
433      /// Constructor
434      CandidateListPivotRule(NetworkSimplex &ns) :
435        _source(ns._source), _target(ns._target),
436        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
437        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
438        _next_arc(0)
439      {
440        // The main parameters of the pivot rule
441        const double LIST_LENGTH_FACTOR = 0.25;
442        const int MIN_LIST_LENGTH = 10;
443        const double MINOR_LIMIT_FACTOR = 0.1;
444        const int MIN_MINOR_LIMIT = 3;
445
446        _list_length = std::max( int(LIST_LENGTH_FACTOR *
447                                     std::sqrt(double(_search_arc_num))),
448                                 MIN_LIST_LENGTH );
449        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
450                                 MIN_MINOR_LIMIT );
451        _curr_length = _minor_count = 0;
452        _candidates.resize(_list_length);
453      }
454
455      /// Find next entering arc
456      bool findEnteringArc() {
457        Cost min, c;
458        int e;
459        if (_curr_length > 0 && _minor_count < _minor_limit) {
460          // Minor iteration: select the best eligible arc from the
461          // current candidate list
462          ++_minor_count;
463          min = 0;
464          for (int i = 0; i < _curr_length; ++i) {
465            e = _candidates[i];
466            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
467            if (c < min) {
468              min = c;
469              _in_arc = e;
470            }
471            else if (c >= 0) {
472              _candidates[i--] = _candidates[--_curr_length];
473            }
474          }
475          if (min < 0) return true;
476        }
477
478        // Major iteration: build a new candidate list
479        min = 0;
480        _curr_length = 0;
481        for (e = _next_arc; e != _search_arc_num; ++e) {
482          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
483          if (c < 0) {
484            _candidates[_curr_length++] = e;
485            if (c < min) {
486              min = c;
487              _in_arc = e;
488            }
489            if (_curr_length == _list_length) goto search_end;
490          }
491        }
492        for (e = 0; e != _next_arc; ++e) {
493          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
494          if (c < 0) {
495            _candidates[_curr_length++] = e;
496            if (c < min) {
497              min = c;
498              _in_arc = e;
499            }
500            if (_curr_length == _list_length) goto search_end;
501          }
502        }
503        if (_curr_length == 0) return false;
504
505      search_end:
506        _minor_count = 1;
507        _next_arc = e;
508        return true;
509      }
510
511    }; //class CandidateListPivotRule
512
513
514    // Implementation of the Altering Candidate List pivot rule
515    class AlteringListPivotRule
516    {
517    private:
518
519      // References to the NetworkSimplex class
520      const IntVector  &_source;
521      const IntVector  &_target;
522      const CostVector &_cost;
523      const CharVector &_state;
524      const CostVector &_pi;
525      int &_in_arc;
526      int _search_arc_num;
527
528      // Pivot rule data
529      int _block_size, _head_length, _curr_length;
530      int _next_arc;
531      IntVector _candidates;
532      CostVector _cand_cost;
533
534      // Functor class to compare arcs during sort of the candidate list
535      class SortFunc
536      {
537      private:
538        const CostVector &_map;
539      public:
540        SortFunc(const CostVector &map) : _map(map) {}
541        bool operator()(int left, int right) {
542          return _map[left] > _map[right];
543        }
544      };
545
546      SortFunc _sort_func;
547
548    public:
549
550      // Constructor
551      AlteringListPivotRule(NetworkSimplex &ns) :
552        _source(ns._source), _target(ns._target),
553        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
554        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
555        _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
556      {
557        // The main parameters of the pivot rule
558        const double BLOCK_SIZE_FACTOR = 1.0;
559        const int MIN_BLOCK_SIZE = 10;
560        const double HEAD_LENGTH_FACTOR = 0.1;
561        const int MIN_HEAD_LENGTH = 3;
562
563        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
564                                    std::sqrt(double(_search_arc_num))),
565                                MIN_BLOCK_SIZE );
566        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
567                                 MIN_HEAD_LENGTH );
568        _candidates.resize(_head_length + _block_size);
569        _curr_length = 0;
570      }
571
572      // Find next entering arc
573      bool findEnteringArc() {
574        // Check the current candidate list
575        int e;
576        Cost c;
577        for (int i = 0; i != _curr_length; ++i) {
578          e = _candidates[i];
579          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
580          if (c < 0) {
581            _cand_cost[e] = c;
582          } else {
583            _candidates[i--] = _candidates[--_curr_length];
584          }
585        }
586
587        // Extend the list
588        int cnt = _block_size;
589        int limit = _head_length;
590
591        for (e = _next_arc; e != _search_arc_num; ++e) {
592          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
593          if (c < 0) {
594            _cand_cost[e] = c;
595            _candidates[_curr_length++] = e;
596          }
597          if (--cnt == 0) {
598            if (_curr_length > limit) goto search_end;
599            limit = 0;
600            cnt = _block_size;
601          }
602        }
603        for (e = 0; e != _next_arc; ++e) {
604          _cand_cost[e] = _state[e] *
605            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
606          if (_cand_cost[e] < 0) {
607            _candidates[_curr_length++] = e;
608          }
609          if (--cnt == 0) {
610            if (_curr_length > limit) goto search_end;
611            limit = 0;
612            cnt = _block_size;
613          }
614        }
615        if (_curr_length == 0) return false;
616
617      search_end:
618
619        // Make heap of the candidate list (approximating a partial sort)
620        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
621                   _sort_func );
622
623        // Pop the first element of the heap
624        _in_arc = _candidates[0];
625        _next_arc = e;
626        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
627                  _sort_func );
628        _curr_length = std::min(_head_length, _curr_length - 1);
629        return true;
630      }
631
632    }; //class AlteringListPivotRule
633
634  public:
635
636    /// \brief Constructor.
637    ///
638    /// The constructor of the class.
639    ///
640    /// \param graph The digraph the algorithm runs on.
641    /// \param arc_mixing Indicate if the arcs will be stored in a
642    /// mixed order in the internal data structure.
643    /// In general, it leads to similar performance as using the original
644    /// arc order, but it makes the algorithm more robust and in special
645    /// cases, even significantly faster. Therefore, it is enabled by default.
646    NetworkSimplex(const GR& graph, bool arc_mixing = true) :
647      _graph(graph), _node_id(graph), _arc_id(graph),
648      _arc_mixing(arc_mixing),
649      MAX(std::numeric_limits<Value>::max()),
650      INF(std::numeric_limits<Value>::has_infinity ?
651          std::numeric_limits<Value>::infinity() : MAX)
652    {
653      // Check the number types
654      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
655        "The flow type of NetworkSimplex must be signed");
656      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
657        "The cost type of NetworkSimplex must be signed");
658
659      // Reset data structures
660      reset();
661    }
662
663    /// \name Parameters
664    /// The parameters of the algorithm can be specified using these
665    /// functions.
666
667    /// @{
668
669    /// \brief Set the lower bounds on the arcs.
670    ///
671    /// This function sets the lower bounds on the arcs.
672    /// If it is not used before calling \ref run(), the lower bounds
673    /// will be set to zero on all arcs.
674    ///
675    /// \param map An arc map storing the lower bounds.
676    /// Its \c Value type must be convertible to the \c Value type
677    /// of the algorithm.
678    ///
679    /// \return <tt>(*this)</tt>
680    template <typename LowerMap>
681    NetworkSimplex& lowerMap(const LowerMap& map) {
682      _have_lower = true;
683      for (ArcIt a(_graph); a != INVALID; ++a) {
684        _lower[_arc_id[a]] = map[a];
685      }
686      return *this;
687    }
688
689    /// \brief Set the upper bounds (capacities) on the arcs.
690    ///
691    /// This function sets the upper bounds (capacities) on the arcs.
692    /// If it is not used before calling \ref run(), the upper bounds
693    /// will be set to \ref INF on all arcs (i.e. the flow value will be
694    /// unbounded from above).
695    ///
696    /// \param map An arc map storing the upper bounds.
697    /// Its \c Value type must be convertible to the \c Value type
698    /// of the algorithm.
699    ///
700    /// \return <tt>(*this)</tt>
701    template<typename UpperMap>
702    NetworkSimplex& upperMap(const UpperMap& map) {
703      for (ArcIt a(_graph); a != INVALID; ++a) {
704        _upper[_arc_id[a]] = map[a];
705      }
706      return *this;
707    }
708
709    /// \brief Set the costs of the arcs.
710    ///
711    /// This function sets the costs of the arcs.
712    /// If it is not used before calling \ref run(), the costs
713    /// will be set to \c 1 on all arcs.
714    ///
715    /// \param map An arc map storing the costs.
716    /// Its \c Value type must be convertible to the \c Cost type
717    /// of the algorithm.
718    ///
719    /// \return <tt>(*this)</tt>
720    template<typename CostMap>
721    NetworkSimplex& costMap(const CostMap& map) {
722      for (ArcIt a(_graph); a != INVALID; ++a) {
723        _cost[_arc_id[a]] = map[a];
724      }
725      return *this;
726    }
727
728    /// \brief Set the supply values of the nodes.
729    ///
730    /// This function sets the supply values of the nodes.
731    /// If neither this function nor \ref stSupply() is used before
732    /// calling \ref run(), the supply of each node will be set to zero.
733    ///
734    /// \param map A node map storing the supply values.
735    /// Its \c Value type must be convertible to the \c Value type
736    /// of the algorithm.
737    ///
738    /// \return <tt>(*this)</tt>
739    ///
740    /// \sa supplyType()
741    template<typename SupplyMap>
742    NetworkSimplex& supplyMap(const SupplyMap& map) {
743      for (NodeIt n(_graph); n != INVALID; ++n) {
744        _supply[_node_id[n]] = map[n];
745      }
746      return *this;
747    }
748
749    /// \brief Set single source and target nodes and a supply value.
750    ///
751    /// This function sets a single source node and a single target node
752    /// and the required flow value.
753    /// If neither this function nor \ref supplyMap() is used before
754    /// calling \ref run(), the supply of each node will be set to zero.
755    ///
756    /// Using this function has the same effect as using \ref supplyMap()
757    /// with a map in which \c k is assigned to \c s, \c -k is
758    /// assigned to \c t and all other nodes have zero supply value.
759    ///
760    /// \param s The source node.
761    /// \param t The target node.
762    /// \param k The required amount of flow from node \c s to node \c t
763    /// (i.e. the supply of \c s and the demand of \c t).
764    ///
765    /// \return <tt>(*this)</tt>
766    NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
767      for (int i = 0; i != _node_num; ++i) {
768        _supply[i] = 0;
769      }
770      _supply[_node_id[s]] =  k;
771      _supply[_node_id[t]] = -k;
772      return *this;
773    }
774
775    /// \brief Set the type of the supply constraints.
776    ///
777    /// This function sets the type of the supply/demand constraints.
778    /// If it is not used before calling \ref run(), the \ref GEQ supply
779    /// type will be used.
780    ///
781    /// For more information, see \ref SupplyType.
782    ///
783    /// \return <tt>(*this)</tt>
784    NetworkSimplex& supplyType(SupplyType supply_type) {
785      _stype = supply_type;
786      return *this;
787    }
788
789    /// @}
790
791    /// \name Execution Control
792    /// The algorithm can be executed using \ref run().
793
794    /// @{
795
796    /// \brief Run the algorithm.
797    ///
798    /// This function runs the algorithm.
799    /// The paramters can be specified using functions \ref lowerMap(),
800    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
801    /// \ref supplyType().
802    /// For example,
803    /// \code
804    ///   NetworkSimplex<ListDigraph> ns(graph);
805    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
806    ///     .supplyMap(sup).run();
807    /// \endcode
808    ///
809    /// This function can be called more than once. All the given parameters
810    /// are kept for the next call, unless \ref resetParams() or \ref reset()
811    /// is used, thus only the modified parameters have to be set again.
812    /// If the underlying digraph was also modified after the construction
813    /// of the class (or the last \ref reset() call), then the \ref reset()
814    /// function must be called.
815    ///
816    /// \param pivot_rule The pivot rule that will be used during the
817    /// algorithm. For more information, see \ref PivotRule.
818    ///
819    /// \return \c INFEASIBLE if no feasible flow exists,
820    /// \n \c OPTIMAL if the problem has optimal solution
821    /// (i.e. it is feasible and bounded), and the algorithm has found
822    /// optimal flow and node potentials (primal and dual solutions),
823    /// \n \c UNBOUNDED if the objective function of the problem is
824    /// unbounded, i.e. there is a directed cycle having negative total
825    /// cost and infinite upper bound.
826    ///
827    /// \see ProblemType, PivotRule
828    /// \see resetParams(), reset()
829    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
830      if (!init()) return INFEASIBLE;
831      return start(pivot_rule);
832    }
833
834    /// \brief Reset all the parameters that have been given before.
835    ///
836    /// This function resets all the paramaters that have been given
837    /// before using functions \ref lowerMap(), \ref upperMap(),
838    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
839    ///
840    /// It is useful for multiple \ref run() calls. Basically, all the given
841    /// parameters are kept for the next \ref run() call, unless
842    /// \ref resetParams() or \ref reset() is used.
843    /// If the underlying digraph was also modified after the construction
844    /// of the class or the last \ref reset() call, then the \ref reset()
845    /// function must be used, otherwise \ref resetParams() is sufficient.
846    ///
847    /// For example,
848    /// \code
849    ///   NetworkSimplex<ListDigraph> ns(graph);
850    ///
851    ///   // First run
852    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
853    ///     .supplyMap(sup).run();
854    ///
855    ///   // Run again with modified cost map (resetParams() is not called,
856    ///   // so only the cost map have to be set again)
857    ///   cost[e] += 100;
858    ///   ns.costMap(cost).run();
859    ///
860    ///   // Run again from scratch using resetParams()
861    ///   // (the lower bounds will be set to zero on all arcs)
862    ///   ns.resetParams();
863    ///   ns.upperMap(capacity).costMap(cost)
864    ///     .supplyMap(sup).run();
865    /// \endcode
866    ///
867    /// \return <tt>(*this)</tt>
868    ///
869    /// \see reset(), run()
870    NetworkSimplex& resetParams() {
871      for (int i = 0; i != _node_num; ++i) {
872        _supply[i] = 0;
873      }
874      for (int i = 0; i != _arc_num; ++i) {
875        _lower[i] = 0;
876        _upper[i] = INF;
877        _cost[i] = 1;
878      }
879      _have_lower = false;
880      _stype = GEQ;
881      return *this;
882    }
883
884    /// \brief Reset the internal data structures and all the parameters
885    /// that have been given before.
886    ///
887    /// This function resets the internal data structures and all the
888    /// paramaters that have been given before using functions \ref lowerMap(),
889    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
890    /// \ref supplyType().
891    ///
892    /// It is useful for multiple \ref run() calls. Basically, all the given
893    /// parameters are kept for the next \ref run() call, unless
894    /// \ref resetParams() or \ref reset() is used.
895    /// If the underlying digraph was also modified after the construction
896    /// of the class or the last \ref reset() call, then the \ref reset()
897    /// function must be used, otherwise \ref resetParams() is sufficient.
898    ///
899    /// See \ref resetParams() for examples.
900    ///
901    /// \return <tt>(*this)</tt>
902    ///
903    /// \see resetParams(), run()
904    NetworkSimplex& reset() {
905      // Resize vectors
906      _node_num = countNodes(_graph);
907      _arc_num = countArcs(_graph);
908      int all_node_num = _node_num + 1;
909      int max_arc_num = _arc_num + 2 * _node_num;
910
911      _source.resize(max_arc_num);
912      _target.resize(max_arc_num);
913
914      _lower.resize(_arc_num);
915      _upper.resize(_arc_num);
916      _cap.resize(max_arc_num);
917      _cost.resize(max_arc_num);
918      _supply.resize(all_node_num);
919      _flow.resize(max_arc_num);
920      _pi.resize(all_node_num);
921
922      _parent.resize(all_node_num);
923      _pred.resize(all_node_num);
924      _pred_dir.resize(all_node_num);
925      _thread.resize(all_node_num);
926      _rev_thread.resize(all_node_num);
927      _succ_num.resize(all_node_num);
928      _last_succ.resize(all_node_num);
929      _state.resize(max_arc_num);
930
931      // Copy the graph
932      int i = 0;
933      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
934        _node_id[n] = i;
935      }
936      if (_arc_mixing) {
937        // Store the arcs in a mixed order
938        const int skip = std::max(_arc_num / _node_num, 3);
939        int i = 0, j = 0;
940        for (ArcIt a(_graph); a != INVALID; ++a) {
941          _arc_id[a] = i;
942          _source[i] = _node_id[_graph.source(a)];
943          _target[i] = _node_id[_graph.target(a)];
944          if ((i += skip) >= _arc_num) i = ++j;
945        }
946      } else {
947        // Store the arcs in the original order
948        int i = 0;
949        for (ArcIt a(_graph); a != INVALID; ++a, ++i) {
950          _arc_id[a] = i;
951          _source[i] = _node_id[_graph.source(a)];
952          _target[i] = _node_id[_graph.target(a)];
953        }
954      }
955
956      // Reset parameters
957      resetParams();
958      return *this;
959    }
960
961    /// @}
962
963    /// \name Query Functions
964    /// The results of the algorithm can be obtained using these
965    /// functions.\n
966    /// The \ref run() function must be called before using them.
967
968    /// @{
969
970    /// \brief Return the total cost of the found flow.
971    ///
972    /// This function returns the total cost of the found flow.
973    /// Its complexity is O(e).
974    ///
975    /// \note The return type of the function can be specified as a
976    /// template parameter. For example,
977    /// \code
978    ///   ns.totalCost<double>();
979    /// \endcode
980    /// It is useful if the total cost cannot be stored in the \c Cost
981    /// type of the algorithm, which is the default return type of the
982    /// function.
983    ///
984    /// \pre \ref run() must be called before using this function.
985    template <typename Number>
986    Number totalCost() const {
987      Number c = 0;
988      for (ArcIt a(_graph); a != INVALID; ++a) {
989        int i = _arc_id[a];
990        c += Number(_flow[i]) * Number(_cost[i]);
991      }
992      return c;
993    }
994
995#ifndef DOXYGEN
996    Cost totalCost() const {
997      return totalCost<Cost>();
998    }
999#endif
1000
1001    /// \brief Return the flow on the given arc.
1002    ///
1003    /// This function returns the flow on the given arc.
1004    ///
1005    /// \pre \ref run() must be called before using this function.
1006    Value flow(const Arc& a) const {
1007      return _flow[_arc_id[a]];
1008    }
1009
1010    /// \brief Copy the flow values (the primal solution) into the
1011    /// given map.
1012    ///
1013    /// This function copies the flow value on each arc into the given
1014    /// map. The \c Value type of the algorithm must be convertible to
1015    /// the \c Value type of the map.
1016    ///
1017    /// \pre \ref run() must be called before using this function.
1018    template <typename FlowMap>
1019    void flowMap(FlowMap &map) const {
1020      for (ArcIt a(_graph); a != INVALID; ++a) {
1021        map.set(a, _flow[_arc_id[a]]);
1022      }
1023    }
1024
1025    /// \brief Return the potential (dual value) of the given node.
1026    ///
1027    /// This function returns the potential (dual value) of the
1028    /// given node.
1029    ///
1030    /// \pre \ref run() must be called before using this function.
1031    Cost potential(const Node& n) const {
1032      return _pi[_node_id[n]];
1033    }
1034
1035    /// \brief Copy the potential values (the dual solution) into the
1036    /// given map.
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.