COIN-OR::LEMON - Graph Library

source: lemon/lemon/network_simplex.h @ 1189:a12cca3ad15a

Last change on this file since 1189:a12cca3ad15a was 1166:d59484d5fc1f, checked in by Alpar Juttner <alpar@…>, 11 years ago

Merge docfix #437

File size: 51.3 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 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(e).
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      // Remove non-zero lower bounds
1071      if (_have_lower) {
1072        for (int i = 0; i != _arc_num; ++i) {
1073          Value c = _lower[i];
1074          if (c >= 0) {
1075            _cap[i] = _upper[i] < MAX ? _upper[i] - c : INF;
1076          } else {
1077            _cap[i] = _upper[i] < MAX + c ? _upper[i] - c : INF;
1078          }
1079          _supply[_source[i]] -= c;
1080          _supply[_target[i]] += c;
1081        }
1082      } else {
1083        for (int i = 0; i != _arc_num; ++i) {
1084          _cap[i] = _upper[i];
1085        }
1086      }
1087
1088      // Initialize artifical cost
1089      Cost ART_COST;
1090      if (std::numeric_limits<Cost>::is_exact) {
1091        ART_COST = std::numeric_limits<Cost>::max() / 2 + 1;
1092      } else {
1093        ART_COST = 0;
1094        for (int i = 0; i != _arc_num; ++i) {
1095          if (_cost[i] > ART_COST) ART_COST = _cost[i];
1096        }
1097        ART_COST = (ART_COST + 1) * _node_num;
1098      }
1099
1100      // Initialize arc maps
1101      for (int i = 0; i != _arc_num; ++i) {
1102        _flow[i] = 0;
1103        _state[i] = STATE_LOWER;
1104      }
1105
1106      // Set data for the artificial root node
1107      _root = _node_num;
1108      _parent[_root] = -1;
1109      _pred[_root] = -1;
1110      _thread[_root] = 0;
1111      _rev_thread[0] = _root;
1112      _succ_num[_root] = _node_num + 1;
1113      _last_succ[_root] = _root - 1;
1114      _supply[_root] = -_sum_supply;
1115      _pi[_root] = 0;
1116
1117      // Add artificial arcs and initialize the spanning tree data structure
1118      if (_sum_supply == 0) {
1119        // EQ supply constraints
1120        _search_arc_num = _arc_num;
1121        _all_arc_num = _arc_num + _node_num;
1122        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1123          _parent[u] = _root;
1124          _pred[u] = e;
1125          _thread[u] = u + 1;
1126          _rev_thread[u + 1] = u;
1127          _succ_num[u] = 1;
1128          _last_succ[u] = u;
1129          _cap[e] = INF;
1130          _state[e] = STATE_TREE;
1131          if (_supply[u] >= 0) {
1132            _pred_dir[u] = DIR_UP;
1133            _pi[u] = 0;
1134            _source[e] = u;
1135            _target[e] = _root;
1136            _flow[e] = _supply[u];
1137            _cost[e] = 0;
1138          } else {
1139            _pred_dir[u] = DIR_DOWN;
1140            _pi[u] = ART_COST;
1141            _source[e] = _root;
1142            _target[e] = u;
1143            _flow[e] = -_supply[u];
1144            _cost[e] = ART_COST;
1145          }
1146        }
1147      }
1148      else if (_sum_supply > 0) {
1149        // LEQ supply constraints
1150        _search_arc_num = _arc_num + _node_num;
1151        int f = _arc_num + _node_num;
1152        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1153          _parent[u] = _root;
1154          _thread[u] = u + 1;
1155          _rev_thread[u + 1] = u;
1156          _succ_num[u] = 1;
1157          _last_succ[u] = u;
1158          if (_supply[u] >= 0) {
1159            _pred_dir[u] = DIR_UP;
1160            _pi[u] = 0;
1161            _pred[u] = e;
1162            _source[e] = u;
1163            _target[e] = _root;
1164            _cap[e] = INF;
1165            _flow[e] = _supply[u];
1166            _cost[e] = 0;
1167            _state[e] = STATE_TREE;
1168          } else {
1169            _pred_dir[u] = DIR_DOWN;
1170            _pi[u] = ART_COST;
1171            _pred[u] = f;
1172            _source[f] = _root;
1173            _target[f] = u;
1174            _cap[f] = INF;
1175            _flow[f] = -_supply[u];
1176            _cost[f] = ART_COST;
1177            _state[f] = STATE_TREE;
1178            _source[e] = u;
1179            _target[e] = _root;
1180            _cap[e] = INF;
1181            _flow[e] = 0;
1182            _cost[e] = 0;
1183            _state[e] = STATE_LOWER;
1184            ++f;
1185          }
1186        }
1187        _all_arc_num = f;
1188      }
1189      else {
1190        // GEQ supply constraints
1191        _search_arc_num = _arc_num + _node_num;
1192        int f = _arc_num + _node_num;
1193        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1194          _parent[u] = _root;
1195          _thread[u] = u + 1;
1196          _rev_thread[u + 1] = u;
1197          _succ_num[u] = 1;
1198          _last_succ[u] = u;
1199          if (_supply[u] <= 0) {
1200            _pred_dir[u] = DIR_DOWN;
1201            _pi[u] = 0;
1202            _pred[u] = e;
1203            _source[e] = _root;
1204            _target[e] = u;
1205            _cap[e] = INF;
1206            _flow[e] = -_supply[u];
1207            _cost[e] = 0;
1208            _state[e] = STATE_TREE;
1209          } else {
1210            _pred_dir[u] = DIR_UP;
1211            _pi[u] = -ART_COST;
1212            _pred[u] = f;
1213            _source[f] = u;
1214            _target[f] = _root;
1215            _cap[f] = INF;
1216            _flow[f] = _supply[u];
1217            _state[f] = STATE_TREE;
1218            _cost[f] = ART_COST;
1219            _source[e] = _root;
1220            _target[e] = u;
1221            _cap[e] = INF;
1222            _flow[e] = 0;
1223            _cost[e] = 0;
1224            _state[e] = STATE_LOWER;
1225            ++f;
1226          }
1227        }
1228        _all_arc_num = f;
1229      }
1230
1231      return true;
1232    }
1233
1234    // Find the join node
1235    void findJoinNode() {
1236      int u = _source[in_arc];
1237      int v = _target[in_arc];
1238      while (u != v) {
1239        if (_succ_num[u] < _succ_num[v]) {
1240          u = _parent[u];
1241        } else {
1242          v = _parent[v];
1243        }
1244      }
1245      join = u;
1246    }
1247
1248    // Find the leaving arc of the cycle and returns true if the
1249    // leaving arc is not the same as the entering arc
1250    bool findLeavingArc() {
1251      // Initialize first and second nodes according to the direction
1252      // of the cycle
1253      int first, second;
1254      if (_state[in_arc] == STATE_LOWER) {
1255        first  = _source[in_arc];
1256        second = _target[in_arc];
1257      } else {
1258        first  = _target[in_arc];
1259        second = _source[in_arc];
1260      }
1261      delta = _cap[in_arc];
1262      int result = 0;
1263      Value c, d;
1264      int e;
1265
1266      // Search the cycle form the first node to the join node
1267      for (int u = first; u != join; u = _parent[u]) {
1268        e = _pred[u];
1269        d = _flow[e];
1270        if (_pred_dir[u] == DIR_DOWN) {
1271          c = _cap[e];
1272          d = c >= MAX ? INF : c - d;
1273        }
1274        if (d < delta) {
1275          delta = d;
1276          u_out = u;
1277          result = 1;
1278        }
1279      }
1280
1281      // Search the cycle form the second node to the join node
1282      for (int u = second; u != join; u = _parent[u]) {
1283        e = _pred[u];
1284        d = _flow[e];
1285        if (_pred_dir[u] == DIR_UP) {
1286          c = _cap[e];
1287          d = c >= MAX ? INF : c - d;
1288        }
1289        if (d <= delta) {
1290          delta = d;
1291          u_out = u;
1292          result = 2;
1293        }
1294      }
1295
1296      if (result == 1) {
1297        u_in = first;
1298        v_in = second;
1299      } else {
1300        u_in = second;
1301        v_in = first;
1302      }
1303      return result != 0;
1304    }
1305
1306    // Change _flow and _state vectors
1307    void changeFlow(bool change) {
1308      // Augment along the cycle
1309      if (delta > 0) {
1310        Value val = _state[in_arc] * delta;
1311        _flow[in_arc] += val;
1312        for (int u = _source[in_arc]; u != join; u = _parent[u]) {
1313          _flow[_pred[u]] -= _pred_dir[u] * val;
1314        }
1315        for (int u = _target[in_arc]; u != join; u = _parent[u]) {
1316          _flow[_pred[u]] += _pred_dir[u] * val;
1317        }
1318      }
1319      // Update the state of the entering and leaving arcs
1320      if (change) {
1321        _state[in_arc] = STATE_TREE;
1322        _state[_pred[u_out]] =
1323          (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
1324      } else {
1325        _state[in_arc] = -_state[in_arc];
1326      }
1327    }
1328
1329    // Update the tree structure
1330    void updateTreeStructure() {
1331      int old_rev_thread = _rev_thread[u_out];
1332      int old_succ_num = _succ_num[u_out];
1333      int old_last_succ = _last_succ[u_out];
1334      v_out = _parent[u_out];
1335
1336      // Check if u_in and u_out coincide
1337      if (u_in == u_out) {
1338        // Update _parent, _pred, _pred_dir
1339        _parent[u_in] = v_in;
1340        _pred[u_in] = in_arc;
1341        _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
1342
1343        // Update _thread and _rev_thread
1344        if (_thread[v_in] != u_out) {
1345          int after = _thread[old_last_succ];
1346          _thread[old_rev_thread] = after;
1347          _rev_thread[after] = old_rev_thread;
1348          after = _thread[v_in];
1349          _thread[v_in] = u_out;
1350          _rev_thread[u_out] = v_in;
1351          _thread[old_last_succ] = after;
1352          _rev_thread[after] = old_last_succ;
1353        }
1354      } else {
1355        // Handle the case when old_rev_thread equals to v_in
1356        // (it also means that join and v_out coincide)
1357        int thread_continue = old_rev_thread == v_in ?
1358          _thread[old_last_succ] : _thread[v_in];
1359
1360        // Update _thread and _parent along the stem nodes (i.e. the nodes
1361        // between u_in and u_out, whose parent have to be changed)
1362        int stem = u_in;              // the current stem node
1363        int par_stem = v_in;          // the new parent of stem
1364        int next_stem;                // the next stem node
1365        int last = _last_succ[u_in];  // the last successor of stem
1366        int before, after = _thread[last];
1367        _thread[v_in] = u_in;
1368        _dirty_revs.clear();
1369        _dirty_revs.push_back(v_in);
1370        while (stem != u_out) {
1371          // Insert the next stem node into the thread list
1372          next_stem = _parent[stem];
1373          _thread[last] = next_stem;
1374          _dirty_revs.push_back(last);
1375
1376          // Remove the subtree of stem from the thread list
1377          before = _rev_thread[stem];
1378          _thread[before] = after;
1379          _rev_thread[after] = before;
1380
1381          // Change the parent node and shift stem nodes
1382          _parent[stem] = par_stem;
1383          par_stem = stem;
1384          stem = next_stem;
1385
1386          // Update last and after
1387          last = _last_succ[stem] == _last_succ[par_stem] ?
1388            _rev_thread[par_stem] : _last_succ[stem];
1389          after = _thread[last];
1390        }
1391        _parent[u_out] = par_stem;
1392        _thread[last] = thread_continue;
1393        _rev_thread[thread_continue] = last;
1394        _last_succ[u_out] = last;
1395
1396        // Remove the subtree of u_out from the thread list except for
1397        // the case when old_rev_thread equals to v_in
1398        if (old_rev_thread != v_in) {
1399          _thread[old_rev_thread] = after;
1400          _rev_thread[after] = old_rev_thread;
1401        }
1402
1403        // Update _rev_thread using the new _thread values
1404        for (int i = 0; i != int(_dirty_revs.size()); ++i) {
1405          int u = _dirty_revs[i];
1406          _rev_thread[_thread[u]] = u;
1407        }
1408
1409        // Update _pred, _pred_dir, _last_succ and _succ_num for the
1410        // stem nodes from u_out to u_in
1411        int tmp_sc = 0, tmp_ls = _last_succ[u_out];
1412        for (int u = u_out, p = _parent[u]; u != u_in; u = p, p = _parent[u]) {
1413          _pred[u] = _pred[p];
1414          _pred_dir[u] = -_pred_dir[p];
1415          tmp_sc += _succ_num[u] - _succ_num[p];
1416          _succ_num[u] = tmp_sc;
1417          _last_succ[p] = tmp_ls;
1418        }
1419        _pred[u_in] = in_arc;
1420        _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
1421        _succ_num[u_in] = old_succ_num;
1422      }
1423
1424      // Update _last_succ from v_in towards the root
1425      int up_limit_out = _last_succ[join] == v_in ? join : -1;
1426      int last_succ_out = _last_succ[u_out];
1427      for (int u = v_in; u != -1 && _last_succ[u] == v_in; u = _parent[u]) {
1428        _last_succ[u] = last_succ_out;
1429      }
1430
1431      // Update _last_succ from v_out towards the root
1432      if (join != old_rev_thread && v_in != old_rev_thread) {
1433        for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1434             u = _parent[u]) {
1435          _last_succ[u] = old_rev_thread;
1436        }
1437      }
1438      else if (last_succ_out != old_last_succ) {
1439        for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1440             u = _parent[u]) {
1441          _last_succ[u] = last_succ_out;
1442        }
1443      }
1444
1445      // Update _succ_num from v_in to join
1446      for (int u = v_in; u != join; u = _parent[u]) {
1447        _succ_num[u] += old_succ_num;
1448      }
1449      // Update _succ_num from v_out to join
1450      for (int u = v_out; u != join; u = _parent[u]) {
1451        _succ_num[u] -= old_succ_num;
1452      }
1453    }
1454
1455    // Update potentials in the subtree that has been moved
1456    void updatePotential() {
1457      Cost sigma = _pi[v_in] - _pi[u_in] -
1458                   _pred_dir[u_in] * _cost[in_arc];
1459      int end = _thread[_last_succ[u_in]];
1460      for (int u = u_in; u != end; u = _thread[u]) {
1461        _pi[u] += sigma;
1462      }
1463    }
1464
1465    // Heuristic initial pivots
1466    bool initialPivots() {
1467      Value curr, total = 0;
1468      std::vector<Node> supply_nodes, demand_nodes;
1469      for (NodeIt u(_graph); u != INVALID; ++u) {
1470        curr = _supply[_node_id[u]];
1471        if (curr > 0) {
1472          total += curr;
1473          supply_nodes.push_back(u);
1474        }
1475        else if (curr < 0) {
1476          demand_nodes.push_back(u);
1477        }
1478      }
1479      if (_sum_supply > 0) total -= _sum_supply;
1480      if (total <= 0) return true;
1481
1482      IntVector arc_vector;
1483      if (_sum_supply >= 0) {
1484        if (supply_nodes.size() == 1 && demand_nodes.size() == 1) {
1485          // Perform a reverse graph search from the sink to the source
1486          typename GR::template NodeMap<bool> reached(_graph, false);
1487          Node s = supply_nodes[0], t = demand_nodes[0];
1488          std::vector<Node> stack;
1489          reached[t] = true;
1490          stack.push_back(t);
1491          while (!stack.empty()) {
1492            Node u, v = stack.back();
1493            stack.pop_back();
1494            if (v == s) break;
1495            for (InArcIt a(_graph, v); a != INVALID; ++a) {
1496              if (reached[u = _graph.source(a)]) continue;
1497              int j = _arc_id[a];
1498              if (_cap[j] >= total) {
1499                arc_vector.push_back(j);
1500                reached[u] = true;
1501                stack.push_back(u);
1502              }
1503            }
1504          }
1505        } else {
1506          // Find the min. cost incomming arc for each demand node
1507          for (int i = 0; i != int(demand_nodes.size()); ++i) {
1508            Node v = demand_nodes[i];
1509            Cost c, min_cost = std::numeric_limits<Cost>::max();
1510            Arc min_arc = INVALID;
1511            for (InArcIt a(_graph, v); a != INVALID; ++a) {
1512              c = _cost[_arc_id[a]];
1513              if (c < min_cost) {
1514                min_cost = c;
1515                min_arc = a;
1516              }
1517            }
1518            if (min_arc != INVALID) {
1519              arc_vector.push_back(_arc_id[min_arc]);
1520            }
1521          }
1522        }
1523      } else {
1524        // Find the min. cost outgoing arc for each supply node
1525        for (int i = 0; i != int(supply_nodes.size()); ++i) {
1526          Node u = supply_nodes[i];
1527          Cost c, min_cost = std::numeric_limits<Cost>::max();
1528          Arc min_arc = INVALID;
1529          for (OutArcIt a(_graph, u); a != INVALID; ++a) {
1530            c = _cost[_arc_id[a]];
1531            if (c < min_cost) {
1532              min_cost = c;
1533              min_arc = a;
1534            }
1535          }
1536          if (min_arc != INVALID) {
1537            arc_vector.push_back(_arc_id[min_arc]);
1538          }
1539        }
1540      }
1541
1542      // Perform heuristic initial pivots
1543      for (int i = 0; i != int(arc_vector.size()); ++i) {
1544        in_arc = arc_vector[i];
1545        if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] -
1546            _pi[_target[in_arc]]) >= 0) continue;
1547        findJoinNode();
1548        bool change = findLeavingArc();
1549        if (delta >= MAX) return false;
1550        changeFlow(change);
1551        if (change) {
1552          updateTreeStructure();
1553          updatePotential();
1554        }
1555      }
1556      return true;
1557    }
1558
1559    // Execute the algorithm
1560    ProblemType start(PivotRule pivot_rule) {
1561      // Select the pivot rule implementation
1562      switch (pivot_rule) {
1563        case FIRST_ELIGIBLE:
1564          return start<FirstEligiblePivotRule>();
1565        case BEST_ELIGIBLE:
1566          return start<BestEligiblePivotRule>();
1567        case BLOCK_SEARCH:
1568          return start<BlockSearchPivotRule>();
1569        case CANDIDATE_LIST:
1570          return start<CandidateListPivotRule>();
1571        case ALTERING_LIST:
1572          return start<AlteringListPivotRule>();
1573      }
1574      return INFEASIBLE; // avoid warning
1575    }
1576
1577    template <typename PivotRuleImpl>
1578    ProblemType start() {
1579      PivotRuleImpl pivot(*this);
1580
1581      // Perform heuristic initial pivots
1582      if (!initialPivots()) return UNBOUNDED;
1583
1584      // Execute the Network Simplex algorithm
1585      while (pivot.findEnteringArc()) {
1586        findJoinNode();
1587        bool change = findLeavingArc();
1588        if (delta >= MAX) return UNBOUNDED;
1589        changeFlow(change);
1590        if (change) {
1591          updateTreeStructure();
1592          updatePotential();
1593        }
1594      }
1595
1596      // Check feasibility
1597      for (int e = _search_arc_num; e != _all_arc_num; ++e) {
1598        if (_flow[e] != 0) return INFEASIBLE;
1599      }
1600
1601      // Transform the solution and the supply map to the original form
1602      if (_have_lower) {
1603        for (int i = 0; i != _arc_num; ++i) {
1604          Value c = _lower[i];
1605          if (c != 0) {
1606            _flow[i] += c;
1607            _supply[_source[i]] += c;
1608            _supply[_target[i]] -= c;
1609          }
1610        }
1611      }
1612
1613      // Shift potentials to meet the requirements of the GEQ/LEQ type
1614      // optimality conditions
1615      if (_sum_supply == 0) {
1616        if (_stype == GEQ) {
1617          Cost max_pot = -std::numeric_limits<Cost>::max();
1618          for (int i = 0; i != _node_num; ++i) {
1619            if (_pi[i] > max_pot) max_pot = _pi[i];
1620          }
1621          if (max_pot > 0) {
1622            for (int i = 0; i != _node_num; ++i)
1623              _pi[i] -= max_pot;
1624          }
1625        } else {
1626          Cost min_pot = std::numeric_limits<Cost>::max();
1627          for (int i = 0; i != _node_num; ++i) {
1628            if (_pi[i] < min_pot) min_pot = _pi[i];
1629          }
1630          if (min_pot < 0) {
1631            for (int i = 0; i != _node_num; ++i)
1632              _pi[i] -= min_pot;
1633          }
1634        }
1635      }
1636
1637      return OPTIMAL;
1638    }
1639
1640  }; //class NetworkSimplex
1641
1642  ///@}
1643
1644} //namespace lemon
1645
1646#endif //LEMON_NETWORK_SIMPLEX_H
Note: See TracBrowser for help on using the repository browser.