COIN-OR::LEMON - Graph Library

source: lemon/lemon/network_simplex.h @ 1269:9eac00ea588f

Last change on this file since 1269:9eac00ea588f was 1254:c5cd8960df74, checked in by Peter Kovacs <kpeter@…>, 6 years ago

Use m instead of e for denoting the number of arcs/edges (#463)

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