COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/network_simplex.h @ 2634:e98bbe64cca4

Last change on this file since 2634:e98bbe64cca4 was 2634:e98bbe64cca4, checked in by Peter Kovacs, 11 years ago

Rework Network Simplex
Use simpler and faster graph implementation instead of SmartGraph?

File size: 37.1 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 = 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 depth vector of the spanning tree structure
595    IntVector _depth;
596    // The parent vector of the spanning tree structure
597    IntVector _parent;
598    // The pred_edge vector of the spanning tree structure
599    IntVector _pred;
600    // The thread vector of the spanning tree structure
601    IntVector _thread;
602    // The forward vector of the spanning tree structure
603    BoolVector _forward;
604    // The state vector
605    IntVector _state;
606    // The root node
607    int _root;
608
609    // The entering edge of the current pivot iteration
610    int _in_edge;
611
612    // Temporary nodes used in the current pivot iteration
613    int join, u_in, v_in, u_out, v_out;
614    int right, first, second, last;
615    int stem, par_stem, new_stem;
616
617    // The maximum augment amount along the found cycle in the current
618    // pivot iteration
619    Capacity delta;
620
621  public:
622
623    /// \brief General constructor (with lower bounds).
624    ///
625    /// General constructor (with lower bounds).
626    ///
627    /// \param graph The directed graph the algorithm runs on.
628    /// \param lower The lower bounds of the edges.
629    /// \param capacity The capacities (upper bounds) of the edges.
630    /// \param cost The cost (length) values of the edges.
631    /// \param supply The supply values of the nodes (signed).
632    NetworkSimplex( const Graph &graph,
633                    const LowerMap &lower,
634                    const CapacityMap &capacity,
635                    const CostMap &cost,
636                    const SupplyMap &supply ) :
637      _orig_graph(graph), _orig_lower(&lower), _orig_cap(capacity),
638      _orig_cost(cost), _orig_supply(&supply),
639      _flow_result(NULL), _potential_result(NULL),
640      _local_flow(false), _local_potential(false),
641      _node_id(graph)
642    {}
643
644    /// \brief General constructor (without lower bounds).
645    ///
646    /// General constructor (without lower bounds).
647    ///
648    /// \param graph The directed graph the algorithm runs on.
649    /// \param capacity The capacities (upper bounds) of the edges.
650    /// \param cost The cost (length) values of the edges.
651    /// \param supply The supply values of the nodes (signed).
652    NetworkSimplex( const Graph &graph,
653                    const CapacityMap &capacity,
654                    const CostMap &cost,
655                    const SupplyMap &supply ) :
656      _orig_graph(graph), _orig_lower(NULL), _orig_cap(capacity),
657      _orig_cost(cost), _orig_supply(&supply),
658      _flow_result(NULL), _potential_result(NULL),
659      _local_flow(false), _local_potential(false),
660      _node_id(graph)
661    {}
662
663    /// \brief Simple constructor (with lower bounds).
664    ///
665    /// Simple constructor (with lower bounds).
666    ///
667    /// \param graph The directed graph the algorithm runs on.
668    /// \param lower The lower bounds of the edges.
669    /// \param capacity The capacities (upper bounds) of the edges.
670    /// \param cost The cost (length) values of the edges.
671    /// \param s The source node.
672    /// \param t The target node.
673    /// \param flow_value The required amount of flow from node \c s
674    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
675    NetworkSimplex( const Graph &graph,
676                    const LowerMap &lower,
677                    const CapacityMap &capacity,
678                    const CostMap &cost,
679                    Node s, Node t,
680                    Capacity flow_value ) :
681      _orig_graph(graph), _orig_lower(&lower), _orig_cap(capacity),
682      _orig_cost(cost), _orig_supply(NULL),
683      _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
684      _flow_result(NULL), _potential_result(NULL),
685      _local_flow(false), _local_potential(false),
686      _node_id(graph)
687    {}
688
689    /// \brief Simple constructor (without lower bounds).
690    ///
691    /// Simple constructor (without lower bounds).
692    ///
693    /// \param graph The directed graph the algorithm runs on.
694    /// \param capacity The capacities (upper bounds) of the edges.
695    /// \param cost The cost (length) values of the edges.
696    /// \param s The source node.
697    /// \param t The target node.
698    /// \param flow_value The required amount of flow from node \c s
699    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
700    NetworkSimplex( const Graph &graph,
701                    const CapacityMap &capacity,
702                    const CostMap &cost,
703                    Node s, Node t,
704                    Capacity flow_value ) :
705      _orig_graph(graph), _orig_lower(NULL), _orig_cap(capacity),
706      _orig_cost(cost), _orig_supply(NULL),
707      _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
708      _flow_result(NULL), _potential_result(NULL),
709      _local_flow(false), _local_potential(false),
710      _node_id(graph)
711    {}
712
713    /// Destructor.
714    ~NetworkSimplex() {
715      if (_local_flow) delete _flow_result;
716      if (_local_potential) delete _potential_result;
717    }
718
719    /// \brief Set the flow map.
720    ///
721    /// Set the flow map.
722    ///
723    /// \return \c (*this)
724    NetworkSimplex& flowMap(FlowMap &map) {
725      if (_local_flow) {
726        delete _flow_result;
727        _local_flow = false;
728      }
729      _flow_result = &map;
730      return *this;
731    }
732
733    /// \brief Set the potential map.
734    ///
735    /// Set the potential map.
736    ///
737    /// \return \c (*this)
738    NetworkSimplex& potentialMap(PotentialMap &map) {
739      if (_local_potential) {
740        delete _potential_result;
741        _local_potential = false;
742      }
743      _potential_result = &map;
744      return *this;
745    }
746
747    /// \name Execution control
748
749    /// @{
750
751    /// \brief Runs the algorithm.
752    ///
753    /// Runs the algorithm.
754    ///
755    /// \param pivot_rule The pivot rule that is used during the
756    /// algorithm.
757    ///
758    /// The available pivot rules:
759    ///
760    /// - FIRST_ELIGIBLE_PIVOT The next eligible edge is selected in
761    /// a wraparound fashion in every iteration
762    /// (\ref FirstEligiblePivotRule).
763    ///
764    /// - BEST_ELIGIBLE_PIVOT The best eligible edge is selected in
765    /// every iteration (\ref BestEligiblePivotRule).
766    ///
767    /// - BLOCK_SEARCH_PIVOT A specified number of edges are examined in
768    /// every iteration in a wraparound fashion and the best eligible
769    /// edge is selected from this block (\ref BlockSearchPivotRule).
770    ///
771    /// - CANDIDATE_LIST_PIVOT In a major iteration a candidate list is
772    /// built from eligible edges in a wraparound fashion and in the
773    /// following minor iterations the best eligible edge is selected
774    /// from this list (\ref CandidateListPivotRule).
775    ///
776    /// - ALTERING_LIST_PIVOT It is a modified version of the
777    /// "Candidate List" pivot rule. It keeps only the several best
778    /// eligible edges from the former candidate list and extends this
779    /// list in every iteration (\ref AlteringListPivotRule).
780    ///
781    /// According to our comprehensive benchmark tests the "Block Search"
782    /// pivot rule proved to be the fastest and the most robust on
783    /// various test inputs. Thus it is the default option.
784    ///
785    /// \return \c true if a feasible flow can be found.
786    bool run(PivotRuleEnum pivot_rule = BLOCK_SEARCH_PIVOT) {
787      return init() && start(pivot_rule);
788    }
789
790    /// @}
791
792    /// \name Query Functions
793    /// The results of the algorithm can be obtained using these
794    /// functions.\n
795    /// \ref lemon::NetworkSimplex::run() "run()" must be called before
796    /// using them.
797
798    /// @{
799
800    /// \brief Return a const reference to the edge map storing the
801    /// found flow.
802    ///
803    /// Return a const reference to the edge map storing the found flow.
804    ///
805    /// \pre \ref run() must be called before using this function.
806    const FlowMap& flowMap() const {
807      return *_flow_result;
808    }
809
810    /// \brief Return a const reference to the node map storing the
811    /// found potentials (the dual solution).
812    ///
813    /// Return a const reference to the node map storing the found
814    /// potentials (the dual solution).
815    ///
816    /// \pre \ref run() must be called before using this function.
817    const PotentialMap& potentialMap() const {
818      return *_potential_result;
819    }
820
821    /// \brief Return the flow on the given edge.
822    ///
823    /// Return the flow on the given edge.
824    ///
825    /// \pre \ref run() must be called before using this function.
826    Capacity flow(const typename Graph::Edge& edge) const {
827      return (*_flow_result)[edge];
828    }
829
830    /// \brief Return the potential of the given node.
831    ///
832    /// Return the potential of the given node.
833    ///
834    /// \pre \ref run() must be called before using this function.
835    Cost potential(const typename Graph::Node& node) const {
836      return (*_potential_result)[node];
837    }
838
839    /// \brief Return the total cost of the found flow.
840    ///
841    /// Return the total cost of the found flow. The complexity of the
842    /// function is \f$ O(e) \f$.
843    ///
844    /// \pre \ref run() must be called before using this function.
845    Cost totalCost() const {
846      Cost c = 0;
847      for (EdgeIt e(_orig_graph); e != INVALID; ++e)
848        c += (*_flow_result)[e] * _orig_cost[e];
849      return c;
850    }
851
852    /// @}
853
854  private:
855
856    // Initialize internal data structures
857    bool init() {
858      // Initialize result maps
859      if (!_flow_result) {
860        _flow_result = new FlowMap(_orig_graph);
861        _local_flow = true;
862      }
863      if (!_potential_result) {
864        _potential_result = new PotentialMap(_orig_graph);
865        _local_potential = true;
866      }
867     
868      // Initialize vectors
869      _node_num = countNodes(_orig_graph);
870      _edge_num = countEdges(_orig_graph);
871      int all_node_num = _node_num + 1;
872      int all_edge_num = _edge_num + _node_num;
873     
874      _edge.resize(_edge_num);
875      _node.reserve(_node_num);
876      _source.resize(all_edge_num);
877      _target.resize(all_edge_num);
878     
879      _cap.resize(all_edge_num);
880      _cost.resize(all_edge_num);
881      _supply.resize(all_node_num);
882      _flow.resize(all_edge_num, 0);
883      _pi.resize(all_node_num, 0);
884     
885      _depth.resize(all_node_num);
886      _parent.resize(all_node_num);
887      _pred.resize(all_node_num);
888      _thread.resize(all_node_num);
889      _forward.resize(all_node_num);
890      _state.resize(all_edge_num, STATE_LOWER);
891     
892      // Initialize node related data
893      bool valid_supply = true;
894      if (_orig_supply) {
895        Supply sum = 0;
896        int i = 0;
897        for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
898          _node.push_back(n);
899          _node_id[n] = i;
900          _supply[i] = (*_orig_supply)[n];
901          sum += _supply[i];
902        }
903        valid_supply = (sum == 0);
904      } else {
905        int i = 0;
906        for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
907          _node.push_back(n);
908          _node_id[n] = i;
909          _supply[i] = 0;
910        }
911        _supply[_node_id[_orig_source]] =  _orig_flow_value;
912        _supply[_node_id[_orig_target]] = -_orig_flow_value;
913      }
914      if (!valid_supply) return false;
915     
916      // Set data for the artificial root node
917      _root = _node_num;
918      _depth[_root] = 0;
919      _parent[_root] = -1;
920      _pred[_root] = -1;
921      _thread[_root] = 0;
922      _supply[_root] = 0;
923      _pi[_root] = 0;
924     
925      // Store the edges in a mixed order
926      int k = std::max(int(sqrt(_edge_num)), 10);
927      int i = 0;
928      for (EdgeIt e(_orig_graph); e != INVALID; ++e) {
929        _edge[i] = e;
930        if ((i += k) >= _edge_num) i = (i % k) + 1;
931      }
932
933      // Initialize edge maps
934      for (int i = 0; i != _edge_num; ++i) {
935        Edge e = _edge[i];
936        _source[i] = _node_id[_orig_graph.source(e)];
937        _target[i] = _node_id[_orig_graph.target(e)];
938        _cost[i] = _orig_cost[e];
939        _cap[i] = _orig_cap[e];
940      }
941
942      // Remove non-zero lower bounds
943      if (_orig_lower) {
944        for (int i = 0; i != _edge_num; ++i) {
945          Capacity c = (*_orig_lower)[_edge[i]];
946          if (c != 0) {
947            _cap[i] -= c;
948            _supply[_source[i]] -= c;
949            _supply[_target[i]] += c;
950          }
951        }
952      }
953
954      // Add artificial edges and initialize the spanning tree data structure
955      Cost max_cost = std::numeric_limits<Cost>::max() / 4;
956      Capacity max_cap = std::numeric_limits<Capacity>::max();
957      for (int u = 0, e = _edge_num; u != _node_num; ++u, ++e) {
958        _thread[u] = u + 1;
959        _depth[u] = 1;
960        _parent[u] = _root;
961        _pred[u] = e;
962        if (_supply[u] >= 0) {
963          _flow[e] = _supply[u];
964          _forward[u] = true;
965          _pi[u] = -max_cost;
966        } else {
967          _flow[e] = -_supply[u];
968          _forward[u] = false;
969          _pi[u] = max_cost;
970        }
971        _cost[e] = max_cost;
972        _cap[e] = max_cap;
973        _state[e] = STATE_TREE;
974      }
975
976      return true;
977    }
978
979    // Find the join node
980    void findJoinNode() {
981      int u = _source[_in_edge];
982      int v = _target[_in_edge];
983      while (_depth[u] > _depth[v]) u = _parent[u];
984      while (_depth[v] > _depth[u]) v = _parent[v];
985      while (u != v) {
986        u = _parent[u];
987        v = _parent[v];
988      }
989      join = u;
990    }
991
992    // Find the leaving edge of the cycle and returns true if the
993    // leaving edge is not the same as the entering edge
994    bool findLeavingEdge() {
995      // Initialize first and second nodes according to the direction
996      // of the cycle
997      if (_state[_in_edge] == STATE_LOWER) {
998        first  = _source[_in_edge];
999        second = _target[_in_edge];
1000      } else {
1001        first  = _target[_in_edge];
1002        second = _source[_in_edge];
1003      }
1004      delta = _cap[_in_edge];
1005      int result = 0;
1006      Capacity d;
1007      int e;
1008
1009      // Search the cycle along the path form the first node to the root
1010      for (int u = first; u != join; u = _parent[u]) {
1011        e = _pred[u];
1012        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
1013        if (d < delta) {
1014          delta = d;
1015          u_out = u;
1016          result = 1;
1017        }
1018      }
1019      // Search the cycle along the path form the second node to the root
1020      for (int u = second; u != join; u = _parent[u]) {
1021        e = _pred[u];
1022        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
1023        if (d <= delta) {
1024          delta = d;
1025          u_out = u;
1026          result = 2;
1027        }
1028      }
1029
1030      if (result == 1) {
1031        u_in = first;
1032        v_in = second;
1033      } else {
1034        u_in = second;
1035        v_in = first;
1036      }
1037      return result != 0;
1038    }
1039
1040    // Change _flow and _state vectors
1041    void changeFlow(bool change) {
1042      // Augment along the cycle
1043      if (delta > 0) {
1044        Capacity val = _state[_in_edge] * delta;
1045        _flow[_in_edge] += val;
1046        for (int u = _source[_in_edge]; u != join; u = _parent[u]) {
1047          _flow[_pred[u]] += _forward[u] ? -val : val;
1048        }
1049        for (int u = _target[_in_edge]; u != join; u = _parent[u]) {
1050          _flow[_pred[u]] += _forward[u] ? val : -val;
1051        }
1052      }
1053      // Update the state of the entering and leaving edges
1054      if (change) {
1055        _state[_in_edge] = STATE_TREE;
1056        _state[_pred[u_out]] =
1057          (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
1058      } else {
1059        _state[_in_edge] = -_state[_in_edge];
1060      }
1061    }
1062
1063    // Update _thread and _parent vectors
1064    void updateThreadParent() {
1065      int u;
1066      v_out = _parent[u_out];
1067
1068      // Handle the case when join and v_out coincide
1069      bool par_first = false;
1070      if (join == v_out) {
1071        for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
1072        if (u == v_in) {
1073          par_first = true;
1074          while (_thread[u] != u_out) u = _thread[u];
1075          first = u;
1076        }
1077      }
1078
1079      // Find the last successor of u_in (u) and the node after it (right)
1080      // according to the thread index
1081      for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
1082      right = _thread[u];
1083      if (_thread[v_in] == u_out) {
1084        for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;
1085        if (last == u_out) last = _thread[last];
1086      }
1087      else last = _thread[v_in];
1088
1089      // Update stem nodes
1090      _thread[v_in] = stem = u_in;
1091      par_stem = v_in;
1092      while (stem != u_out) {
1093        _thread[u] = new_stem = _parent[stem];
1094
1095        // Find the node just before the stem node (u) according to
1096        // the original thread index
1097        for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
1098        _thread[u] = right;
1099
1100        // Change the parent node of stem and shift stem and par_stem nodes
1101        _parent[stem] = par_stem;
1102        par_stem = stem;
1103        stem = new_stem;
1104
1105        // Find the last successor of stem (u) and the node after it (right)
1106        // according to the thread index
1107        for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
1108        right = _thread[u];
1109      }
1110      _parent[u_out] = par_stem;
1111      _thread[u] = last;
1112
1113      if (join == v_out && par_first) {
1114        if (first != v_in) _thread[first] = right;
1115      } else {
1116        for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
1117        _thread[u] = right;
1118      }
1119    }
1120
1121    // Update _pred and _forward vectors
1122    void updatePredEdge() {
1123      int u = u_out, v;
1124      while (u != u_in) {
1125        v = _parent[u];
1126        _pred[u] = _pred[v];
1127        _forward[u] = !_forward[v];
1128        u = v;
1129      }
1130      _pred[u_in] = _in_edge;
1131      _forward[u_in] = (u_in == _source[_in_edge]);
1132    }
1133
1134    // Update _depth and _potential vectors
1135    void updateDepthPotential() {
1136      _depth[u_in] = _depth[v_in] + 1;
1137      Cost sigma = _forward[u_in] ?
1138        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
1139        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
1140      _pi[u_in] += sigma;
1141      for(int u = _thread[u_in]; _parent[u] != -1; u = _thread[u]) {
1142        _depth[u] = _depth[_parent[u]] + 1;
1143        if (_depth[u] <= _depth[u_in]) break;
1144        _pi[u] += sigma;
1145      }
1146    }
1147
1148    // Execute the algorithm
1149    bool start(PivotRuleEnum pivot_rule) {
1150      // Select the pivot rule implementation
1151      switch (pivot_rule) {
1152        case FIRST_ELIGIBLE_PIVOT:
1153          return start<FirstEligiblePivotRule>();
1154        case BEST_ELIGIBLE_PIVOT:
1155          return start<BestEligiblePivotRule>();
1156        case BLOCK_SEARCH_PIVOT:
1157          return start<BlockSearchPivotRule>();
1158        case CANDIDATE_LIST_PIVOT:
1159          return start<CandidateListPivotRule>();
1160        case ALTERING_LIST_PIVOT:
1161          return start<AlteringListPivotRule>();
1162      }
1163      return false;
1164    }
1165
1166    template<class PivotRuleImplementation>
1167    bool start() {
1168      PivotRuleImplementation pivot(*this);
1169
1170      // Execute the network simplex algorithm
1171      while (pivot.findEnteringEdge()) {
1172        findJoinNode();
1173        bool change = findLeavingEdge();
1174        changeFlow(change);
1175        if (change) {
1176          updateThreadParent();
1177          updatePredEdge();
1178          updateDepthPotential();
1179        }
1180      }
1181
1182      // Check if the flow amount equals zero on all the artificial edges
1183      for (int e = _edge_num; e != _edge_num + _node_num; ++e) {
1184        if (_flow[e] > 0) return false;
1185      }
1186
1187      // Copy flow values to _flow_result
1188      if (_orig_lower) {
1189        for (int i = 0; i != _edge_num; ++i) {
1190          Edge e = _edge[i];
1191          (*_flow_result)[e] = (*_orig_lower)[e] + _flow[i];
1192        }
1193      } else {
1194        for (int i = 0; i != _edge_num; ++i) {
1195          (*_flow_result)[_edge[i]] = _flow[i];
1196        }
1197      }
1198      // Copy potential values to _potential_result
1199      for (int i = 0; i != _node_num; ++i) {
1200        (*_potential_result)[_node[i]] = _pi[i];
1201      }
1202
1203      return true;
1204    }
1205
1206  }; //class NetworkSimplex
1207
1208  ///@}
1209
1210} //namespace lemon
1211
1212#endif //LEMON_NETWORK_SIMPLEX_H
Note: See TracBrowser for help on using the repository browser.