COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/network_simplex.h @ 2635:23570e368afa

Last change on this file since 2635:23570e368afa was 2635:23570e368afa, checked in by Peter Kovacs, 10 years ago

XTI data structure for NS - backport hg commit [8c3112a66878]

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