COIN-OR::LEMON - Graph Library

source: lemon/lemon/network_simplex.h @ 655:6ac5d9ae1d3d

Last change on this file since 655:6ac5d9ae1d3d was 655:6ac5d9ae1d3d, checked in by Peter Kovacs <kpeter@…>, 11 years ago

Support real types + numerical stability fix in NS (#254)

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