COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/network_simplex.h

Last change on this file was 2638:61bf01404ede, checked in by Peter Kovacs, 10 years ago

Various improvements for NS pivot rules

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