COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/network_simplex.h @ 605:5232721b3f14

Last change on this file since 605:5232721b3f14 was 605:5232721b3f14, checked in by Peter Kovacs <kpeter@…>, 11 years ago

Rework the interface of NetworkSimplex? (#234)

The parameters of the problem can be set with separate functions
instead of different constructors.

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