COIN-OR::LEMON - Graph Library

source: lemon/lemon/network_simplex.h @ 648:e8349c6f12ca

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

Port NetworkSimplex? from SVN -r3520 (#234)

File size: 36.1 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/math.h>
32
33namespace lemon {
34
35  /// \addtogroup min_cost_flow
36  /// @{
37
38  /// \brief Implementation of the primal network simplex algorithm
39  /// for finding a \ref min_cost_flow "minimum cost flow".
40  ///
41  /// \ref NetworkSimplex implements the primal network simplex algorithm
42  /// for finding a \ref min_cost_flow "minimum cost flow".
43  ///
44  /// \tparam Digraph The digraph type the algorithm runs on.
45  /// \tparam LowerMap The type of the lower bound map.
46  /// \tparam CapacityMap The type of the capacity (upper bound) map.
47  /// \tparam CostMap The type of the cost (length) map.
48  /// \tparam SupplyMap The type of the supply map.
49  ///
50  /// \warning
51  /// - Arc capacities and costs should be \e non-negative \e integers.
52  /// - Supply values should be \e signed \e integers.
53  /// - The value types of the maps should be convertible to each other.
54  /// - \c CostMap::Value must be signed type.
55  ///
56  /// \note \ref NetworkSimplex provides five different pivot rule
57  /// implementations that significantly affect the efficiency of the
58  /// algorithm.
59  /// By default "Block Search" pivot rule is used, which proved to be
60  /// by far the most efficient according to our benchmark tests.
61  /// However another pivot rule can be selected using \ref run()
62  /// function with the proper parameter.
63#ifdef DOXYGEN
64  template < typename Digraph,
65             typename LowerMap,
66             typename CapacityMap,
67             typename CostMap,
68             typename SupplyMap >
69
70#else
71  template < typename Digraph,
72             typename LowerMap = typename Digraph::template ArcMap<int>,
73             typename CapacityMap = typename Digraph::template ArcMap<int>,
74             typename CostMap = typename Digraph::template ArcMap<int>,
75             typename SupplyMap = typename Digraph::template NodeMap<int> >
76#endif
77  class NetworkSimplex
78  {
79    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
80
81    typedef typename CapacityMap::Value Capacity;
82    typedef typename CostMap::Value Cost;
83    typedef typename SupplyMap::Value Supply;
84
85    typedef std::vector<Arc> ArcVector;
86    typedef std::vector<Node> NodeVector;
87    typedef std::vector<int> IntVector;
88    typedef std::vector<bool> BoolVector;
89    typedef std::vector<Capacity> CapacityVector;
90    typedef std::vector<Cost> CostVector;
91    typedef std::vector<Supply> SupplyVector;
92
93  public:
94
95    /// The type of the flow map
96    typedef typename Digraph::template ArcMap<Capacity> FlowMap;
97    /// The type of the potential map
98    typedef typename Digraph::template NodeMap<Cost> PotentialMap;
99
100  public:
101
102    /// Enum type for selecting the pivot rule used by \ref run()
103    enum PivotRuleEnum {
104      FIRST_ELIGIBLE_PIVOT,
105      BEST_ELIGIBLE_PIVOT,
106      BLOCK_SEARCH_PIVOT,
107      CANDIDATE_LIST_PIVOT,
108      ALTERING_LIST_PIVOT
109    };
110
111  private:
112
113    // State constants for arcs
114    enum ArcStateEnum {
115      STATE_UPPER = -1,
116      STATE_TREE  =  0,
117      STATE_LOWER =  1
118    };
119
120  private:
121
122    // References for the original data
123    const Digraph &_orig_graph;
124    const LowerMap *_orig_lower;
125    const CapacityMap &_orig_cap;
126    const CostMap &_orig_cost;
127    const SupplyMap *_orig_supply;
128    Node _orig_source;
129    Node _orig_target;
130    Capacity _orig_flow_value;
131
132    // Result maps
133    FlowMap *_flow_result;
134    PotentialMap *_potential_result;
135    bool _local_flow;
136    bool _local_potential;
137
138    // Data structures for storing the graph
139    ArcVector _arc;
140    NodeVector _node;
141    IntNodeMap _node_id;
142    IntVector _source;
143    IntVector _target;
144
145    // The number of nodes and arcs in the original graph
146    int _node_num;
147    int _arc_num;
148
149    // Node and arc maps
150    CapacityVector _cap;
151    CostVector _cost;
152    CostVector _supply;
153    CapacityVector _flow;
154    CostVector _pi;
155
156    // Node and arc maps for the spanning tree structure
157    IntVector _depth;
158    IntVector _parent;
159    IntVector _pred;
160    IntVector _thread;
161    BoolVector _forward;
162    IntVector _state;
163
164    // The root node
165    int _root;
166
167    // The entering arc in the current pivot iteration
168    int _in_arc;
169
170    // Temporary data used in the current pivot iteration
171    int join, u_in, v_in, u_out, v_out;
172    int right, first, second, last;
173    int stem, par_stem, new_stem;
174    Capacity delta;
175
176  private:
177
178    /// \brief Implementation of the "First Eligible" pivot rule for the
179    /// \ref NetworkSimplex "network simplex" algorithm.
180    ///
181    /// This class implements the "First Eligible" pivot rule
182    /// for the \ref NetworkSimplex "network simplex" algorithm.
183    ///
184    /// For more information see \ref NetworkSimplex::run().
185    class FirstEligiblePivotRule
186    {
187    private:
188
189      // References to the NetworkSimplex class
190      const ArcVector &_arc;
191      const IntVector  &_source;
192      const IntVector  &_target;
193      const CostVector &_cost;
194      const IntVector  &_state;
195      const CostVector &_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        _arc(ns._arc), _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        Cost 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    /// \brief Implementation of the "Best Eligible" pivot rule for the
237    /// \ref NetworkSimplex "network simplex" algorithm.
238    ///
239    /// This class implements the "Best Eligible" pivot rule
240    /// for the \ref NetworkSimplex "network simplex" algorithm.
241    ///
242    /// For more information see \ref NetworkSimplex::run().
243    class BestEligiblePivotRule
244    {
245    private:
246
247      // References to the NetworkSimplex class
248      const ArcVector &_arc;
249      const IntVector  &_source;
250      const IntVector  &_target;
251      const CostVector &_cost;
252      const IntVector  &_state;
253      const CostVector &_pi;
254      int &_in_arc;
255      int _arc_num;
256
257    public:
258
259      /// Constructor
260      BestEligiblePivotRule(NetworkSimplex &ns) :
261        _arc(ns._arc), _source(ns._source), _target(ns._target),
262        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
263        _in_arc(ns._in_arc), _arc_num(ns._arc_num)
264      {}
265
266      /// Find next entering arc
267      bool findEnteringArc() {
268        Cost c, min = 0;
269        for (int e = 0; e < _arc_num; ++e) {
270          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
271          if (c < min) {
272            min = c;
273            _in_arc = e;
274          }
275        }
276        return min < 0;
277      }
278
279    }; //class BestEligiblePivotRule
280
281
282    /// \brief Implementation of the "Block Search" pivot rule for the
283    /// \ref NetworkSimplex "network simplex" algorithm.
284    ///
285    /// This class implements the "Block Search" pivot rule
286    /// for the \ref NetworkSimplex "network simplex" algorithm.
287    ///
288    /// For more information see \ref NetworkSimplex::run().
289    class BlockSearchPivotRule
290    {
291    private:
292
293      // References to the NetworkSimplex class
294      const ArcVector &_arc;
295      const IntVector  &_source;
296      const IntVector  &_target;
297      const CostVector &_cost;
298      const IntVector  &_state;
299      const CostVector &_pi;
300      int &_in_arc;
301      int _arc_num;
302
303      // Pivot rule data
304      int _block_size;
305      int _next_arc;
306
307    public:
308
309      /// Constructor
310      BlockSearchPivotRule(NetworkSimplex &ns) :
311        _arc(ns._arc), _source(ns._source), _target(ns._target),
312        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
313        _in_arc(ns._in_arc), _arc_num(ns._arc_num), _next_arc(0)
314      {
315        // The main parameters of the pivot rule
316        const double BLOCK_SIZE_FACTOR = 2.0;
317        const int MIN_BLOCK_SIZE = 10;
318
319        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
320                                MIN_BLOCK_SIZE );
321      }
322
323      /// Find next entering arc
324      bool findEnteringArc() {
325        Cost c, min = 0;
326        int cnt = _block_size;
327        int e, min_arc = _next_arc;
328        for (e = _next_arc; e < _arc_num; ++e) {
329          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
330          if (c < min) {
331            min = c;
332            min_arc = e;
333          }
334          if (--cnt == 0) {
335            if (min < 0) break;
336            cnt = _block_size;
337          }
338        }
339        if (min == 0 || cnt > 0) {
340          for (e = 0; e < _next_arc; ++e) {
341            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
342            if (c < min) {
343              min = c;
344              min_arc = e;
345            }
346            if (--cnt == 0) {
347              if (min < 0) break;
348              cnt = _block_size;
349            }
350          }
351        }
352        if (min >= 0) return false;
353        _in_arc = min_arc;
354        _next_arc = e;
355        return true;
356      }
357
358    }; //class BlockSearchPivotRule
359
360
361    /// \brief Implementation of the "Candidate List" pivot rule for the
362    /// \ref NetworkSimplex "network simplex" algorithm.
363    ///
364    /// This class implements the "Candidate List" pivot rule
365    /// for the \ref NetworkSimplex "network simplex" algorithm.
366    ///
367    /// For more information see \ref NetworkSimplex::run().
368    class CandidateListPivotRule
369    {
370    private:
371
372      // References to the NetworkSimplex class
373      const ArcVector &_arc;
374      const IntVector  &_source;
375      const IntVector  &_target;
376      const CostVector &_cost;
377      const IntVector  &_state;
378      const CostVector &_pi;
379      int &_in_arc;
380      int _arc_num;
381
382      // Pivot rule data
383      IntVector _candidates;
384      int _list_length, _minor_limit;
385      int _curr_length, _minor_count;
386      int _next_arc;
387
388    public:
389
390      /// Constructor
391      CandidateListPivotRule(NetworkSimplex &ns) :
392        _arc(ns._arc), _source(ns._source), _target(ns._target),
393        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
394        _in_arc(ns._in_arc), _arc_num(ns._arc_num), _next_arc(0)
395      {
396        // The main parameters of the pivot rule
397        const double LIST_LENGTH_FACTOR = 1.0;
398        const int MIN_LIST_LENGTH = 10;
399        const double MINOR_LIMIT_FACTOR = 0.1;
400        const int MIN_MINOR_LIMIT = 3;
401
402        _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_arc_num)),
403                                 MIN_LIST_LENGTH );
404        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
405                                 MIN_MINOR_LIMIT );
406        _curr_length = _minor_count = 0;
407        _candidates.resize(_list_length);
408      }
409
410      /// Find next entering arc
411      bool findEnteringArc() {
412        Cost min, c;
413        int e, min_arc = _next_arc;
414        if (_curr_length > 0 && _minor_count < _minor_limit) {
415          // Minor iteration: select the best eligible arc from the
416          // current candidate list
417          ++_minor_count;
418          min = 0;
419          for (int i = 0; i < _curr_length; ++i) {
420            e = _candidates[i];
421            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
422            if (c < min) {
423              min = c;
424              min_arc = e;
425            }
426            if (c >= 0) {
427              _candidates[i--] = _candidates[--_curr_length];
428            }
429          }
430          if (min < 0) {
431            _in_arc = min_arc;
432            return true;
433          }
434        }
435
436        // Major iteration: build a new candidate list
437        min = 0;
438        _curr_length = 0;
439        for (e = _next_arc; e < _arc_num; ++e) {
440          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
441          if (c < 0) {
442            _candidates[_curr_length++] = e;
443            if (c < min) {
444              min = c;
445              min_arc = e;
446            }
447            if (_curr_length == _list_length) break;
448          }
449        }
450        if (_curr_length < _list_length) {
451          for (e = 0; e < _next_arc; ++e) {
452            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
453            if (c < 0) {
454              _candidates[_curr_length++] = e;
455              if (c < min) {
456                min = c;
457                min_arc = e;
458              }
459              if (_curr_length == _list_length) break;
460            }
461          }
462        }
463        if (_curr_length == 0) return false;
464        _minor_count = 1;
465        _in_arc = min_arc;
466        _next_arc = e;
467        return true;
468      }
469
470    }; //class CandidateListPivotRule
471
472
473    /// \brief Implementation of the "Altering Candidate List" pivot rule
474    /// for the \ref NetworkSimplex "network simplex" algorithm.
475    ///
476    /// This class implements the "Altering Candidate List" pivot rule
477    /// for the \ref NetworkSimplex "network simplex" algorithm.
478    ///
479    /// For more information see \ref NetworkSimplex::run().
480    class AlteringListPivotRule
481    {
482    private:
483
484      // References to the NetworkSimplex class
485      const ArcVector &_arc;
486      const IntVector  &_source;
487      const IntVector  &_target;
488      const CostVector &_cost;
489      const IntVector  &_state;
490      const CostVector &_pi;
491      int &_in_arc;
492      int _arc_num;
493
494      // Pivot rule data
495      int _block_size, _head_length, _curr_length;
496      int _next_arc;
497      IntVector _candidates;
498      CostVector _cand_cost;
499
500      // Functor class to compare arcs during sort of the candidate list
501      class SortFunc
502      {
503      private:
504        const CostVector &_map;
505      public:
506        SortFunc(const CostVector &map) : _map(map) {}
507        bool operator()(int left, int right) {
508          return _map[left] > _map[right];
509        }
510      };
511
512      SortFunc _sort_func;
513
514    public:
515
516      /// Constructor
517      AlteringListPivotRule(NetworkSimplex &ns) :
518        _arc(ns._arc), _source(ns._source), _target(ns._target),
519        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
520        _in_arc(ns._in_arc), _arc_num(ns._arc_num),
521        _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
522      {
523        // The main parameters of the pivot rule
524        const double BLOCK_SIZE_FACTOR = 1.5;
525        const int MIN_BLOCK_SIZE = 10;
526        const double HEAD_LENGTH_FACTOR = 0.1;
527        const int MIN_HEAD_LENGTH = 3;
528
529        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
530                                MIN_BLOCK_SIZE );
531        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
532                                 MIN_HEAD_LENGTH );
533        _candidates.resize(_head_length + _block_size);
534        _curr_length = 0;
535      }
536
537      /// Find next entering arc
538      bool findEnteringArc() {
539        // Check the current candidate list
540        int e;
541        for (int i = 0; i < _curr_length; ++i) {
542          e = _candidates[i];
543          _cand_cost[e] = _state[e] *
544            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
545          if (_cand_cost[e] >= 0) {
546            _candidates[i--] = _candidates[--_curr_length];
547          }
548        }
549
550        // Extend the list
551        int cnt = _block_size;
552        int last_edge = 0;
553        int limit = _head_length;
554
555        for (int e = _next_arc; e < _arc_num; ++e) {
556          _cand_cost[e] = _state[e] *
557            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
558          if (_cand_cost[e] < 0) {
559            _candidates[_curr_length++] = e;
560            last_edge = e;
561          }
562          if (--cnt == 0) {
563            if (_curr_length > limit) break;
564            limit = 0;
565            cnt = _block_size;
566          }
567        }
568        if (_curr_length <= limit) {
569          for (int e = 0; e < _next_arc; ++e) {
570            _cand_cost[e] = _state[e] *
571              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
572            if (_cand_cost[e] < 0) {
573              _candidates[_curr_length++] = e;
574              last_edge = e;
575            }
576            if (--cnt == 0) {
577              if (_curr_length > limit) break;
578              limit = 0;
579              cnt = _block_size;
580            }
581          }
582        }
583        if (_curr_length == 0) return false;
584        _next_arc = last_edge + 1;
585
586        // Make heap of the candidate list (approximating a partial sort)
587        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
588                   _sort_func );
589
590        // Pop the first element of the heap
591        _in_arc = _candidates[0];
592        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
593                  _sort_func );
594        _curr_length = std::min(_head_length, _curr_length - 1);
595        return true;
596      }
597
598    }; //class AlteringListPivotRule
599
600  public:
601
602    /// \brief General constructor (with lower bounds).
603    ///
604    /// General constructor (with lower bounds).
605    ///
606    /// \param digraph The digraph the algorithm runs on.
607    /// \param lower The lower bounds of the arcs.
608    /// \param capacity The capacities (upper bounds) of the arcs.
609    /// \param cost The cost (length) values of the arcs.
610    /// \param supply The supply values of the nodes (signed).
611    NetworkSimplex( const Digraph &digraph,
612                    const LowerMap &lower,
613                    const CapacityMap &capacity,
614                    const CostMap &cost,
615                    const SupplyMap &supply ) :
616      _orig_graph(digraph), _orig_lower(&lower), _orig_cap(capacity),
617      _orig_cost(cost), _orig_supply(&supply),
618      _flow_result(NULL), _potential_result(NULL),
619      _local_flow(false), _local_potential(false),
620      _node_id(digraph)
621    {}
622
623    /// \brief General constructor (without lower bounds).
624    ///
625    /// General constructor (without lower bounds).
626    ///
627    /// \param digraph The digraph the algorithm runs on.
628    /// \param capacity The capacities (upper bounds) of the arcs.
629    /// \param cost The cost (length) values of the arcs.
630    /// \param supply The supply values of the nodes (signed).
631    NetworkSimplex( const Digraph &digraph,
632                    const CapacityMap &capacity,
633                    const CostMap &cost,
634                    const SupplyMap &supply ) :
635      _orig_graph(digraph), _orig_lower(NULL), _orig_cap(capacity),
636      _orig_cost(cost), _orig_supply(&supply),
637      _flow_result(NULL), _potential_result(NULL),
638      _local_flow(false), _local_potential(false),
639      _node_id(digraph)
640    {}
641
642    /// \brief Simple constructor (with lower bounds).
643    ///
644    /// Simple constructor (with lower bounds).
645    ///
646    /// \param digraph The digraph the algorithm runs on.
647    /// \param lower The lower bounds of the arcs.
648    /// \param capacity The capacities (upper bounds) of the arcs.
649    /// \param cost The cost (length) values of the arcs.
650    /// \param s The source node.
651    /// \param t The target node.
652    /// \param flow_value The required amount of flow from node \c s
653    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
654    NetworkSimplex( const Digraph &digraph,
655                    const LowerMap &lower,
656                    const CapacityMap &capacity,
657                    const CostMap &cost,
658                    Node s, Node t,
659                    Capacity flow_value ) :
660      _orig_graph(digraph), _orig_lower(&lower), _orig_cap(capacity),
661      _orig_cost(cost), _orig_supply(NULL),
662      _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
663      _flow_result(NULL), _potential_result(NULL),
664      _local_flow(false), _local_potential(false),
665      _node_id(digraph)
666    {}
667
668    /// \brief Simple constructor (without lower bounds).
669    ///
670    /// Simple constructor (without lower bounds).
671    ///
672    /// \param digraph The digraph the algorithm runs on.
673    /// \param capacity The capacities (upper bounds) of the arcs.
674    /// \param cost The cost (length) values of the arcs.
675    /// \param s The source node.
676    /// \param t The target node.
677    /// \param flow_value The required amount of flow from node \c s
678    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
679    NetworkSimplex( const Digraph &digraph,
680                    const CapacityMap &capacity,
681                    const CostMap &cost,
682                    Node s, Node t,
683                    Capacity flow_value ) :
684      _orig_graph(digraph), _orig_lower(NULL), _orig_cap(capacity),
685      _orig_cost(cost), _orig_supply(NULL),
686      _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
687      _flow_result(NULL), _potential_result(NULL),
688      _local_flow(false), _local_potential(false),
689      _node_id(digraph)
690    {}
691
692    /// Destructor.
693    ~NetworkSimplex() {
694      if (_local_flow) delete _flow_result;
695      if (_local_potential) delete _potential_result;
696    }
697
698    /// \brief Set the flow map.
699    ///
700    /// This function sets the flow map.
701    ///
702    /// \return <tt>(*this)</tt>
703    NetworkSimplex& flowMap(FlowMap &map) {
704      if (_local_flow) {
705        delete _flow_result;
706        _local_flow = false;
707      }
708      _flow_result = &map;
709      return *this;
710    }
711
712    /// \brief Set the potential map.
713    ///
714    /// This function sets the potential map.
715    ///
716    /// \return <tt>(*this)</tt>
717    NetworkSimplex& potentialMap(PotentialMap &map) {
718      if (_local_potential) {
719        delete _potential_result;
720        _local_potential = false;
721      }
722      _potential_result = &map;
723      return *this;
724    }
725
726    /// \name Execution control
727    /// The algorithm can be executed using the
728    /// \ref lemon::NetworkSimplex::run() "run()" function.
729    /// @{
730
731    /// \brief Run the algorithm.
732    ///
733    /// This function runs the algorithm.
734    ///
735    /// \param pivot_rule The pivot rule that is used during the
736    /// algorithm.
737    ///
738    /// The available pivot rules:
739    ///
740    /// - FIRST_ELIGIBLE_PIVOT The next eligible arc is selected in
741    /// a wraparound fashion in every iteration
742    /// (\ref FirstEligiblePivotRule).
743    ///
744    /// - BEST_ELIGIBLE_PIVOT The best eligible arc is selected in
745    /// every iteration (\ref BestEligiblePivotRule).
746    ///
747    /// - BLOCK_SEARCH_PIVOT A specified number of arcs are examined in
748    /// every iteration in a wraparound fashion and the best eligible
749    /// arc is selected from this block (\ref BlockSearchPivotRule).
750    ///
751    /// - CANDIDATE_LIST_PIVOT In a major iteration a candidate list is
752    /// built from eligible arcs in a wraparound fashion and in the
753    /// following minor iterations the best eligible arc is selected
754    /// from this list (\ref CandidateListPivotRule).
755    ///
756    /// - ALTERING_LIST_PIVOT It is a modified version of the
757    /// "Candidate List" pivot rule. It keeps only the several best
758    /// eligible arcs from the former candidate list and extends this
759    /// list in every iteration (\ref AlteringListPivotRule).
760    ///
761    /// According to our comprehensive benchmark tests the "Block Search"
762    /// pivot rule proved to be the fastest and the most robust on
763    /// various test inputs. Thus it is the default option.
764    ///
765    /// \return \c true if a feasible flow can be found.
766    bool run(PivotRuleEnum pivot_rule = BLOCK_SEARCH_PIVOT) {
767      return init() && start(pivot_rule);
768    }
769
770    /// @}
771
772    /// \name Query Functions
773    /// The results of the algorithm can be obtained using these
774    /// functions.\n
775    /// \ref lemon::NetworkSimplex::run() "run()" must be called before
776    /// using them.
777    /// @{
778
779    /// \brief Return a const reference to the flow map.
780    ///
781    /// This function returns a const reference to an arc map storing
782    /// the found flow.
783    ///
784    /// \pre \ref run() must be called before using this function.
785    const FlowMap& flowMap() const {
786      return *_flow_result;
787    }
788
789    /// \brief Return a const reference to the potential map
790    /// (the dual solution).
791    ///
792    /// This function returns a const reference to a node map storing
793    /// the found potentials (the dual solution).
794    ///
795    /// \pre \ref run() must be called before using this function.
796    const PotentialMap& potentialMap() const {
797      return *_potential_result;
798    }
799
800    /// \brief Return the flow on the given arc.
801    ///
802    /// This function returns the flow on the given arc.
803    ///
804    /// \pre \ref run() must be called before using this function.
805    Capacity flow(const Arc& arc) const {
806      return (*_flow_result)[arc];
807    }
808
809    /// \brief Return the potential of the given node.
810    ///
811    /// This function returns the potential of the given node.
812    ///
813    /// \pre \ref run() must be called before using this function.
814    Cost potential(const Node& node) const {
815      return (*_potential_result)[node];
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    /// \pre \ref run() must be called before using this function.
824    Cost totalCost() const {
825      Cost c = 0;
826      for (ArcIt e(_orig_graph); e != INVALID; ++e)
827        c += (*_flow_result)[e] * _orig_cost[e];
828      return c;
829    }
830
831    /// @}
832
833  private:
834
835    // Initialize internal data structures
836    bool init() {
837      // Initialize result maps
838      if (!_flow_result) {
839        _flow_result = new FlowMap(_orig_graph);
840        _local_flow = true;
841      }
842      if (!_potential_result) {
843        _potential_result = new PotentialMap(_orig_graph);
844        _local_potential = true;
845      }
846
847      // Initialize vectors
848      _node_num = countNodes(_orig_graph);
849      _arc_num = countArcs(_orig_graph);
850      int all_node_num = _node_num + 1;
851      int all_edge_num = _arc_num + _node_num;
852
853      _arc.resize(_arc_num);
854      _node.reserve(_node_num);
855      _source.resize(all_edge_num);
856      _target.resize(all_edge_num);
857
858      _cap.resize(all_edge_num);
859      _cost.resize(all_edge_num);
860      _supply.resize(all_node_num);
861      _flow.resize(all_edge_num, 0);
862      _pi.resize(all_node_num, 0);
863
864      _depth.resize(all_node_num);
865      _parent.resize(all_node_num);
866      _pred.resize(all_node_num);
867      _thread.resize(all_node_num);
868      _forward.resize(all_node_num);
869      _state.resize(all_edge_num, STATE_LOWER);
870
871      // Initialize node related data
872      bool valid_supply = true;
873      if (_orig_supply) {
874        Supply sum = 0;
875        int i = 0;
876        for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
877          _node.push_back(n);
878          _node_id[n] = i;
879          _supply[i] = (*_orig_supply)[n];
880          sum += _supply[i];
881        }
882        valid_supply = (sum == 0);
883      } else {
884        int i = 0;
885        for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
886          _node.push_back(n);
887          _node_id[n] = i;
888          _supply[i] = 0;
889        }
890        _supply[_node_id[_orig_source]] =  _orig_flow_value;
891        _supply[_node_id[_orig_target]] = -_orig_flow_value;
892      }
893      if (!valid_supply) return false;
894
895      // Set data for the artificial root node
896      _root = _node_num;
897      _depth[_root] = 0;
898      _parent[_root] = -1;
899      _pred[_root] = -1;
900      _thread[_root] = 0;
901      _supply[_root] = 0;
902      _pi[_root] = 0;
903
904      // Store the arcs in a mixed order
905      int k = std::max(int(sqrt(_arc_num)), 10);
906      int i = 0;
907      for (ArcIt e(_orig_graph); e != INVALID; ++e) {
908        _arc[i] = e;
909        if ((i += k) >= _arc_num) i = (i % k) + 1;
910      }
911
912      // Initialize arc maps
913      for (int i = 0; i != _arc_num; ++i) {
914        Arc e = _arc[i];
915        _source[i] = _node_id[_orig_graph.source(e)];
916        _target[i] = _node_id[_orig_graph.target(e)];
917        _cost[i] = _orig_cost[e];
918        _cap[i] = _orig_cap[e];
919      }
920
921      // Remove non-zero lower bounds
922      if (_orig_lower) {
923        for (int i = 0; i != _arc_num; ++i) {
924          Capacity c = (*_orig_lower)[_arc[i]];
925          if (c != 0) {
926            _cap[i] -= c;
927            _supply[_source[i]] -= c;
928            _supply[_target[i]] += c;
929          }
930        }
931      }
932
933      // Add artificial arcs and initialize the spanning tree data structure
934      Cost max_cost = std::numeric_limits<Cost>::max() / 4;
935      Capacity max_cap = std::numeric_limits<Capacity>::max();
936      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
937        _thread[u] = u + 1;
938        _depth[u] = 1;
939        _parent[u] = _root;
940        _pred[u] = e;
941        if (_supply[u] >= 0) {
942          _flow[e] = _supply[u];
943          _forward[u] = true;
944          _pi[u] = -max_cost;
945        } else {
946          _flow[e] = -_supply[u];
947          _forward[u] = false;
948          _pi[u] = max_cost;
949        }
950        _cost[e] = max_cost;
951        _cap[e] = max_cap;
952        _state[e] = STATE_TREE;
953      }
954
955      return true;
956    }
957
958    // Find the join node
959    void findJoinNode() {
960      int u = _source[_in_arc];
961      int v = _target[_in_arc];
962      while (_depth[u] > _depth[v]) u = _parent[u];
963      while (_depth[v] > _depth[u]) v = _parent[v];
964      while (u != v) {
965        u = _parent[u];
966        v = _parent[v];
967      }
968      join = u;
969    }
970
971    // Find the leaving arc of the cycle and returns true if the
972    // leaving arc is not the same as the entering arc
973    bool findLeavingArc() {
974      // Initialize first and second nodes according to the direction
975      // of the cycle
976      if (_state[_in_arc] == STATE_LOWER) {
977        first  = _source[_in_arc];
978        second = _target[_in_arc];
979      } else {
980        first  = _target[_in_arc];
981        second = _source[_in_arc];
982      }
983      delta = _cap[_in_arc];
984      int result = 0;
985      Capacity d;
986      int e;
987
988      // Search the cycle along the path form the first node to the root
989      for (int u = first; u != join; u = _parent[u]) {
990        e = _pred[u];
991        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
992        if (d < delta) {
993          delta = d;
994          u_out = u;
995          result = 1;
996        }
997      }
998      // Search the cycle along the path form the second node to the root
999      for (int u = second; u != join; u = _parent[u]) {
1000        e = _pred[u];
1001        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
1002        if (d <= delta) {
1003          delta = d;
1004          u_out = u;
1005          result = 2;
1006        }
1007      }
1008
1009      if (result == 1) {
1010        u_in = first;
1011        v_in = second;
1012      } else {
1013        u_in = second;
1014        v_in = first;
1015      }
1016      return result != 0;
1017    }
1018
1019    // Change _flow and _state vectors
1020    void changeFlow(bool change) {
1021      // Augment along the cycle
1022      if (delta > 0) {
1023        Capacity val = _state[_in_arc] * delta;
1024        _flow[_in_arc] += val;
1025        for (int u = _source[_in_arc]; u != join; u = _parent[u]) {
1026          _flow[_pred[u]] += _forward[u] ? -val : val;
1027        }
1028        for (int u = _target[_in_arc]; u != join; u = _parent[u]) {
1029          _flow[_pred[u]] += _forward[u] ? val : -val;
1030        }
1031      }
1032      // Update the state of the entering and leaving arcs
1033      if (change) {
1034        _state[_in_arc] = STATE_TREE;
1035        _state[_pred[u_out]] =
1036          (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
1037      } else {
1038        _state[_in_arc] = -_state[_in_arc];
1039      }
1040    }
1041
1042    // Update _thread and _parent vectors
1043    void updateThreadParent() {
1044      int u;
1045      v_out = _parent[u_out];
1046
1047      // Handle the case when join and v_out coincide
1048      bool par_first = false;
1049      if (join == v_out) {
1050        for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
1051        if (u == v_in) {
1052          par_first = true;
1053          while (_thread[u] != u_out) u = _thread[u];
1054          first = u;
1055        }
1056      }
1057
1058      // Find the last successor of u_in (u) and the node after it (right)
1059      // according to the thread index
1060      for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
1061      right = _thread[u];
1062      if (_thread[v_in] == u_out) {
1063        for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;
1064        if (last == u_out) last = _thread[last];
1065      }
1066      else last = _thread[v_in];
1067
1068      // Update stem nodes
1069      _thread[v_in] = stem = u_in;
1070      par_stem = v_in;
1071      while (stem != u_out) {
1072        _thread[u] = new_stem = _parent[stem];
1073
1074        // Find the node just before the stem node (u) according to
1075        // the original thread index
1076        for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
1077        _thread[u] = right;
1078
1079        // Change the parent node of stem and shift stem and par_stem nodes
1080        _parent[stem] = par_stem;
1081        par_stem = stem;
1082        stem = new_stem;
1083
1084        // Find the last successor of stem (u) and the node after it (right)
1085        // according to the thread index
1086        for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
1087        right = _thread[u];
1088      }
1089      _parent[u_out] = par_stem;
1090      _thread[u] = last;
1091
1092      if (join == v_out && par_first) {
1093        if (first != v_in) _thread[first] = right;
1094      } else {
1095        for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
1096        _thread[u] = right;
1097      }
1098    }
1099
1100    // Update _pred and _forward vectors
1101    void updatePredArc() {
1102      int u = u_out, v;
1103      while (u != u_in) {
1104        v = _parent[u];
1105        _pred[u] = _pred[v];
1106        _forward[u] = !_forward[v];
1107        u = v;
1108      }
1109      _pred[u_in] = _in_arc;
1110      _forward[u_in] = (u_in == _source[_in_arc]);
1111    }
1112
1113    // Update _depth and _potential vectors
1114    void updateDepthPotential() {
1115      _depth[u_in] = _depth[v_in] + 1;
1116      Cost sigma = _forward[u_in] ?
1117        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
1118        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
1119      _pi[u_in] += sigma;
1120      for(int u = _thread[u_in]; _parent[u] != -1; u = _thread[u]) {
1121        _depth[u] = _depth[_parent[u]] + 1;
1122        if (_depth[u] <= _depth[u_in]) break;
1123        _pi[u] += sigma;
1124      }
1125    }
1126
1127    // Execute the algorithm
1128    bool start(PivotRuleEnum pivot_rule) {
1129      // Select the pivot rule implementation
1130      switch (pivot_rule) {
1131        case FIRST_ELIGIBLE_PIVOT:
1132          return start<FirstEligiblePivotRule>();
1133        case BEST_ELIGIBLE_PIVOT:
1134          return start<BestEligiblePivotRule>();
1135        case BLOCK_SEARCH_PIVOT:
1136          return start<BlockSearchPivotRule>();
1137        case CANDIDATE_LIST_PIVOT:
1138          return start<CandidateListPivotRule>();
1139        case ALTERING_LIST_PIVOT:
1140          return start<AlteringListPivotRule>();
1141      }
1142      return false;
1143    }
1144
1145    template<class PivotRuleImplementation>
1146    bool start() {
1147      PivotRuleImplementation pivot(*this);
1148
1149      // Execute the network simplex algorithm
1150      while (pivot.findEnteringArc()) {
1151        findJoinNode();
1152        bool change = findLeavingArc();
1153        changeFlow(change);
1154        if (change) {
1155          updateThreadParent();
1156          updatePredArc();
1157          updateDepthPotential();
1158        }
1159      }
1160
1161      // Check if the flow amount equals zero on all the artificial arcs
1162      for (int e = _arc_num; e != _arc_num + _node_num; ++e) {
1163        if (_flow[e] > 0) return false;
1164      }
1165
1166      // Copy flow values to _flow_result
1167      if (_orig_lower) {
1168        for (int i = 0; i != _arc_num; ++i) {
1169          Arc e = _arc[i];
1170          (*_flow_result)[e] = (*_orig_lower)[e] + _flow[i];
1171        }
1172      } else {
1173        for (int i = 0; i != _arc_num; ++i) {
1174          (*_flow_result)[_arc[i]] = _flow[i];
1175        }
1176      }
1177      // Copy potential values to _potential_result
1178      for (int i = 0; i != _node_num; ++i) {
1179        (*_potential_result)[_node[i]] = _pi[i];
1180      }
1181
1182      return true;
1183    }
1184
1185  }; //class NetworkSimplex
1186
1187  ///@}
1188
1189} //namespace lemon
1190
1191#endif //LEMON_NETWORK_SIMPLEX_H
Note: See TracBrowser for help on using the repository browser.