COIN-OR::LEMON - Graph Library

source: lemon/lemon/network_simplex.h @ 653:c7d160f73d52

Last change on this file since 653:c7d160f73d52 was 653:c7d160f73d52, checked in by Peter Kovacs <kpeter@…>, 10 years ago

Support multiple run() calls in NetworkSimplex? (#234)

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