COIN-OR::LEMON - Graph Library

source: lemon/lemon/network_simplex.h @ 651:8c3112a66878

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

Use XTI implementation instead of ATI in NetworkSimplex? (#234)

XTI (eXtended Threaded Index) is an imporved version of the widely
known ATI (Augmented Threaded Index) method for storing and updating
the spanning tree structure in Network Simplex algorithms.

In the ATI data structure three indices are stored for each node:
predecessor, thread and depth. In the XTI data structure depth is
replaced by the number of successors and the last successor
(according to the thread index).

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