COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/network_simplex.h @ 2619:30fb4d68b0e8

Last change on this file since 2619:30fb4d68b0e8 was 2619:30fb4d68b0e8, checked in by Peter Kovacs, 11 years ago

Improve network simplex algorithm

  • Remove "Limited Search" and "Combined" pivot rules.
  • Add a new pivot rule "Altering Candidate List".
  • Make the edge selection faster in every pivot rule.
  • Set the default rule to "Block Search".
  • Doc improvements.

The algorithm became about 15-35 percent faster on various input files.
"Block Search" pivot rule proved to be by far the fastest on all inputs.

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