COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/network_simplex.h @ 2579:691ce54544c5

Last change on this file since 2579:691ce54544c5 was 2579:691ce54544c5, checked in by Peter Kovacs, 13 years ago

Bug fixes in min cost flow files.
Use enum type instead of static constants in NetworkSimplex? to avoid
linker errors.

File size: 33.4 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
35namespace lemon {
36
37  /// \addtogroup min_cost_flow
38  /// @{
39
40  /// \brief Implementation of the network simplex algorithm for
41  /// finding a minimum cost flow.
42  ///
43  /// \ref NetworkSimplex implements the network simplex algorithm for
44  /// finding a minimum cost flow.
45  ///
46  /// \tparam Graph The directed graph type the algorithm runs on.
47  /// \tparam LowerMap The type of the lower bound map.
48  /// \tparam CapacityMap The type of the capacity (upper bound) map.
49  /// \tparam CostMap The type of the cost (length) map.
50  /// \tparam SupplyMap The type of the supply map.
51  ///
52  /// \warning
53  /// - Edge capacities and costs should be \e non-negative \e integers.
54  /// - Supply values should be \e signed \e integers.
55  /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.
56  /// - \c CapacityMap::Value and \c SupplyMap::Value must be
57  ///   convertible to each other.
58  /// - All value types must be convertible to \c CostMap::Value, which
59  ///   must be signed type.
60  ///
61  /// \note \ref NetworkSimplex provides six different pivot rule
62  /// implementations that significantly affect the efficiency of the
63  /// algorithm.
64  /// By default a combined pivot rule is used, which is the fastest
65  /// implementation according to our benchmark tests.
66  /// Another pivot rule can be selected using \ref run() function
67  /// with the proper parameter.
68  ///
69  /// \author Peter Kovacs
70
71  template < typename Graph,
72             typename LowerMap = typename Graph::template EdgeMap<int>,
73             typename CapacityMap = typename Graph::template EdgeMap<int>,
74             typename CostMap = typename Graph::template EdgeMap<int>,
75             typename SupplyMap = typename Graph::template NodeMap<int> >
76  class NetworkSimplex
77  {
78    typedef typename CapacityMap::Value Capacity;
79    typedef typename CostMap::Value Cost;
80    typedef typename SupplyMap::Value Supply;
81
82    typedef SmartGraph SGraph;
83    GRAPH_TYPEDEFS(typename SGraph);
84
85    typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap;
86    typedef typename SGraph::template EdgeMap<Cost> SCostMap;
87    typedef typename SGraph::template NodeMap<Supply> SSupplyMap;
88    typedef typename SGraph::template NodeMap<Cost> SPotentialMap;
89
90    typedef typename SGraph::template NodeMap<int> IntNodeMap;
91    typedef typename SGraph::template NodeMap<bool> BoolNodeMap;
92    typedef typename SGraph::template NodeMap<Node> NodeNodeMap;
93    typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap;
94    typedef typename SGraph::template EdgeMap<int> IntEdgeMap;
95
96    typedef typename Graph::template NodeMap<Node> NodeRefMap;
97    typedef typename Graph::template EdgeMap<Edge> EdgeRefMap;
98
99  public:
100
101    /// The type of the flow map.
102    typedef typename Graph::template EdgeMap<Capacity> FlowMap;
103    /// The type of the potential map.
104    typedef typename Graph::template NodeMap<Cost> PotentialMap;
105
106  public:
107
108    /// Enum type to select the pivot rule used by \ref run().
109    enum PivotRuleEnum {
110      FIRST_ELIGIBLE_PIVOT,
111      BEST_ELIGIBLE_PIVOT,
112      BLOCK_SEARCH_PIVOT,
113      LIMITED_SEARCH_PIVOT,
114      CANDIDATE_LIST_PIVOT,
115      COMBINED_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    class FirstEligiblePivotRule
155    {
156    private:
157
158      NetworkSimplex &_ns;
159      EdgeIt _next_edge;
160
161    public:
162
163      /// Constructor.
164      FirstEligiblePivotRule(NetworkSimplex &ns) :
165        _ns(ns), _next_edge(ns._graph) {}
166
167      /// Finds the next entering edge.
168      bool findEnteringEdge() {
169        for (EdgeIt e = _next_edge; e != INVALID; ++e) {
170          if (_ns._state[e] * _ns._red_cost[e] < 0) {
171            _ns._in_edge = e;
172            _next_edge = ++e;
173            return true;
174          }
175        }
176        for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
177          if (_ns._state[e] * _ns._red_cost[e] < 0) {
178            _ns._in_edge = e;
179            _next_edge = ++e;
180            return true;
181          }
182        }
183        return false;
184      }
185    }; //class FirstEligiblePivotRule
186
187    /// \brief Implementation of the "Best Eligible" pivot rule for the
188    /// \ref NetworkSimplex "network simplex" algorithm.
189    ///
190    /// This class implements the "Best Eligible" pivot rule
191    /// for the \ref NetworkSimplex "network simplex" algorithm.
192    class BestEligiblePivotRule
193    {
194    private:
195
196      NetworkSimplex &_ns;
197
198    public:
199
200      /// Constructor.
201      BestEligiblePivotRule(NetworkSimplex &ns) : _ns(ns) {}
202
203      /// Finds the next entering edge.
204      bool findEnteringEdge() {
205        Cost min = 0;
206        for (EdgeIt e(_ns._graph); e != INVALID; ++e) {
207          if (_ns._state[e] * _ns._red_cost[e] < min) {
208            min = _ns._state[e] * _ns._red_cost[e];
209            _ns._in_edge = e;
210          }
211        }
212        return min < 0;
213      }
214    }; //class BestEligiblePivotRule
215
216    /// \brief Implementation of the "Block Search" pivot rule for the
217    /// \ref NetworkSimplex "network simplex" algorithm.
218    ///
219    /// This class implements the "Block Search" pivot rule
220    /// for the \ref NetworkSimplex "network simplex" algorithm.
221    class BlockSearchPivotRule
222    {
223    private:
224
225      NetworkSimplex &_ns;
226      EdgeIt _next_edge, _min_edge;
227      int _block_size;
228
229      static const int MIN_BLOCK_SIZE = 10;
230
231    public:
232
233      /// Constructor.
234      BlockSearchPivotRule(NetworkSimplex &ns) :
235        _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
236      {
237        _block_size = 2 * int(sqrt(countEdges(_ns._graph)));
238        if (_block_size < MIN_BLOCK_SIZE) _block_size = MIN_BLOCK_SIZE;
239      }
240
241      /// Finds the next entering edge.
242      bool findEnteringEdge() {
243        Cost curr, min = 0;
244        int cnt = 0;
245        for (EdgeIt e = _next_edge; e != INVALID; ++e) {
246          if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
247            min = curr;
248            _min_edge = e;
249          }
250          if (++cnt == _block_size) {
251            if (min < 0) break;
252            cnt = 0;
253          }
254        }
255        if (min == 0) {
256          for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
257            if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
258              min = curr;
259              _min_edge = e;
260            }
261            if (++cnt == _block_size) {
262              if (min < 0) break;
263              cnt = 0;
264            }
265          }
266        }
267        _ns._in_edge = _min_edge;
268        _next_edge = ++_min_edge;
269        return min < 0;
270      }
271    }; //class BlockSearchPivotRule
272
273    /// \brief Implementation of the "Limited Search" pivot rule for the
274    /// \ref NetworkSimplex "network simplex" algorithm.
275    ///
276    /// This class implements the "Limited Search" pivot rule
277    /// for the \ref NetworkSimplex "network simplex" algorithm.
278    class LimitedSearchPivotRule
279    {
280    private:
281
282      NetworkSimplex &_ns;
283      EdgeIt _next_edge, _min_edge;
284      int _sample_size;
285
286      static const int MIN_SAMPLE_SIZE = 10;
287      static const double SAMPLE_SIZE_FACTOR = 0.0015;
288
289    public:
290
291      /// Constructor.
292      LimitedSearchPivotRule(NetworkSimplex &ns) :
293        _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
294      {
295        _sample_size = int(SAMPLE_SIZE_FACTOR * countEdges(_ns._graph));
296        if (_sample_size < MIN_SAMPLE_SIZE)
297          _sample_size = MIN_SAMPLE_SIZE;
298      }
299
300      /// Finds the next entering edge.
301      bool findEnteringEdge() {
302        Cost curr, min = 0;
303        int cnt = 0;
304        for (EdgeIt e = _next_edge; e != INVALID; ++e) {
305          if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
306            min = curr;
307            _min_edge = e;
308          }
309          if (curr < 0 && ++cnt == _sample_size) break;
310        }
311        if (min == 0) {
312          for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
313            if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
314              min = curr;
315              _min_edge = e;
316            }
317            if (curr < 0 && ++cnt == _sample_size) break;
318          }
319        }
320        _ns._in_edge = _min_edge;
321        _next_edge = ++_min_edge;
322        return min < 0;
323      }
324    }; //class LimitedSearchPivotRule
325
326    /// \brief Implementation of the "Candidate List" pivot rule for the
327    /// \ref NetworkSimplex "network simplex" algorithm.
328    ///
329    /// This class implements the "Candidate List" pivot rule
330    /// for the \ref NetworkSimplex "network simplex" algorithm.
331    class CandidateListPivotRule
332    {
333    private:
334
335      NetworkSimplex &_ns;
336
337      // The list of candidate edges.
338      std::vector<Edge> _candidates;
339      // The maximum length of the edge list.
340      int _list_length;
341      // The maximum number of minor iterations between two major
342      // itarations.
343      int _minor_limit;
344
345      int _minor_count;
346      EdgeIt _next_edge;
347
348      static const double LIST_LENGTH_FACTOR = 0.002;
349      static const double MINOR_LIMIT_FACTOR = 0.1;
350      static const int MIN_LIST_LENGTH = 10;
351      static const int MIN_MINOR_LIMIT = 2;
352
353    public:
354
355      /// Constructor.
356      CandidateListPivotRule(NetworkSimplex &ns) :
357        _ns(ns), _next_edge(ns._graph)
358      {
359        int edge_num = countEdges(_ns._graph);
360        _minor_count = 0;
361        _list_length = int(edge_num * LIST_LENGTH_FACTOR);
362        if (_list_length < MIN_LIST_LENGTH)
363          _list_length = MIN_LIST_LENGTH;
364        _minor_limit = int(_list_length * MINOR_LIMIT_FACTOR);
365        if (_minor_limit < MIN_MINOR_LIMIT)
366          _minor_limit = MIN_MINOR_LIMIT;
367      }
368
369      /// Finds the next entering edge.
370      bool findEnteringEdge() {
371        Cost min, curr;
372        if (_minor_count < _minor_limit && _candidates.size() > 0) {
373          // Minor iteration
374          ++_minor_count;
375          Edge e;
376          min = 0;
377          for (int i = 0; i < int(_candidates.size()); ++i) {
378            e = _candidates[i];
379            if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
380              min = curr;
381              _ns._in_edge = e;
382            }
383          }
384          if (min < 0) return true;
385        }
386
387        // Major iteration
388        _candidates.clear();
389        EdgeIt e = _next_edge;
390        min = 0;
391        for ( ; e != INVALID; ++e) {
392          if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
393            _candidates.push_back(e);
394            if (curr < min) {
395              min = curr;
396              _ns._in_edge = e;
397            }
398            if (int(_candidates.size()) == _list_length) break;
399          }
400        }
401        if (int(_candidates.size()) < _list_length) {
402          for (e = EdgeIt(_ns._graph); e != _next_edge; ++e) {
403            if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
404              _candidates.push_back(e);
405              if (curr < min) {
406                min = curr;
407                _ns._in_edge = e;
408              }
409              if (int(_candidates.size()) == _list_length) break;
410            }
411          }
412        }
413        if (_candidates.size() == 0) return false;
414        _minor_count = 1;
415        _next_edge = ++e;
416        return true;
417      }
418    }; //class CandidateListPivotRule
419
420  private:
421
422    // State constants for edges
423    enum EdgeStateEnum {
424      STATE_UPPER = -1,
425      STATE_TREE  =  0,
426      STATE_LOWER =  1
427    };
428
429    // Constant for the combined pivot rule.
430    static const int COMBINED_PIVOT_MAX_DEG = 5;
431
432  private:
433
434    // The directed graph the algorithm runs on
435    SGraph _graph;
436    // The original graph
437    const Graph &_graph_ref;
438    // The original lower bound map
439    const LowerMap *_lower;
440    // The capacity map
441    SCapacityMap _capacity;
442    // The cost map
443    SCostMap _cost;
444    // The supply map
445    SSupplyMap _supply;
446    bool _valid_supply;
447
448    // Edge map of the current flow
449    SCapacityMap _flow;
450    // Node map of the current potentials
451    SPotentialMap _potential;
452
453    // The depth node map of the spanning tree structure
454    IntNodeMap _depth;
455    // The parent node map of the spanning tree structure
456    NodeNodeMap _parent;
457    // The pred_edge node map of the spanning tree structure
458    EdgeNodeMap _pred_edge;
459    // The thread node map of the spanning tree structure
460    NodeNodeMap _thread;
461    // The forward node map of the spanning tree structure
462    BoolNodeMap _forward;
463    // The state edge map
464    IntEdgeMap _state;
465    // The root node of the starting spanning tree
466    Node _root;
467
468    // The reduced cost map
469    ReducedCostMap _red_cost;
470
471    // Members for handling the original graph
472    FlowMap _flow_result;
473    PotentialMap _potential_result;
474    NodeRefMap _node_ref;
475    EdgeRefMap _edge_ref;
476
477    // The entering edge of the current pivot iteration.
478    Edge _in_edge;
479
480    // Temporary nodes used in the current pivot iteration.
481    Node join, u_in, v_in, u_out, v_out;
482    Node right, first, second, last;
483    Node stem, par_stem, new_stem;
484    // The maximum augment amount along the found cycle in the current
485    // pivot iteration.
486    Capacity delta;
487
488  public :
489
490    /// \brief General constructor of the class (with lower bounds).
491    ///
492    /// General constructor of the class (with lower bounds).
493    ///
494    /// \param graph The directed graph the algorithm runs on.
495    /// \param lower The lower bounds of the edges.
496    /// \param capacity The capacities (upper bounds) of the edges.
497    /// \param cost The cost (length) values of the edges.
498    /// \param supply The supply values of the nodes (signed).
499    NetworkSimplex( const Graph &graph,
500                    const LowerMap &lower,
501                    const CapacityMap &capacity,
502                    const CostMap &cost,
503                    const SupplyMap &supply ) :
504      _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
505      _cost(_graph), _supply(_graph), _flow(_graph),
506      _potential(_graph), _depth(_graph), _parent(_graph),
507      _pred_edge(_graph), _thread(_graph), _forward(_graph),
508      _state(_graph), _red_cost(_graph, _cost, _potential),
509      _flow_result(graph), _potential_result(graph),
510      _node_ref(graph), _edge_ref(graph)
511    {
512      // Checking the sum of supply values
513      Supply sum = 0;
514      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
515        sum += supply[n];
516      if (!(_valid_supply = sum == 0)) return;
517
518      // Copying _graph_ref to _graph
519      _graph.reserveNode(countNodes(_graph_ref) + 1);
520      _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
521      copyGraph(_graph, _graph_ref)
522        .edgeMap(_cost, cost)
523        .nodeRef(_node_ref)
524        .edgeRef(_edge_ref)
525        .run();
526
527      // Removing non-zero lower bounds
528      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
529        _capacity[_edge_ref[e]] = capacity[e] - lower[e];
530      }
531      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {
532        Supply s = supply[n];
533        for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)
534          s += lower[e];
535        for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)
536          s -= lower[e];
537        _supply[_node_ref[n]] = s;
538      }
539    }
540
541    /// \brief General constructor of the class (without lower bounds).
542    ///
543    /// General constructor of the class (without lower bounds).
544    ///
545    /// \param graph The directed graph the algorithm runs on.
546    /// \param capacity The capacities (upper bounds) of the edges.
547    /// \param cost The cost (length) values of the edges.
548    /// \param supply The supply values of the nodes (signed).
549    NetworkSimplex( const Graph &graph,
550                    const CapacityMap &capacity,
551                    const CostMap &cost,
552                    const SupplyMap &supply ) :
553      _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
554      _cost(_graph), _supply(_graph), _flow(_graph),
555      _potential(_graph), _depth(_graph), _parent(_graph),
556      _pred_edge(_graph), _thread(_graph), _forward(_graph),
557      _state(_graph), _red_cost(_graph, _cost, _potential),
558      _flow_result(graph), _potential_result(graph),
559      _node_ref(graph), _edge_ref(graph)
560    {
561      // Checking the sum of supply values
562      Supply sum = 0;
563      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
564        sum += supply[n];
565      if (!(_valid_supply = sum == 0)) return;
566
567      // Copying _graph_ref to graph
568      copyGraph(_graph, _graph_ref)
569        .edgeMap(_capacity, capacity)
570        .edgeMap(_cost, cost)
571        .nodeMap(_supply, supply)
572        .nodeRef(_node_ref)
573        .edgeRef(_edge_ref)
574        .run();
575    }
576
577    /// \brief Simple constructor of the class (with lower bounds).
578    ///
579    /// Simple constructor of the class (with lower bounds).
580    ///
581    /// \param graph The directed graph the algorithm runs on.
582    /// \param lower The lower bounds of the edges.
583    /// \param capacity The capacities (upper bounds) of the edges.
584    /// \param cost The cost (length) values of the edges.
585    /// \param s The source node.
586    /// \param t The target node.
587    /// \param flow_value The required amount of flow from node \c s
588    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
589    NetworkSimplex( const Graph &graph,
590                    const LowerMap &lower,
591                    const CapacityMap &capacity,
592                    const CostMap &cost,
593                    typename Graph::Node s,
594                    typename Graph::Node t,
595                    typename SupplyMap::Value flow_value ) :
596      _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
597      _cost(_graph), _supply(_graph), _flow(_graph),
598      _potential(_graph), _depth(_graph), _parent(_graph),
599      _pred_edge(_graph), _thread(_graph), _forward(_graph),
600      _state(_graph), _red_cost(_graph, _cost, _potential),
601      _flow_result(graph), _potential_result(graph),
602      _node_ref(graph), _edge_ref(graph)
603    {
604      // Copying _graph_ref to graph
605      copyGraph(_graph, _graph_ref)
606        .edgeMap(_cost, cost)
607        .nodeRef(_node_ref)
608        .edgeRef(_edge_ref)
609        .run();
610
611      // Removing non-zero lower bounds
612      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
613        _capacity[_edge_ref[e]] = capacity[e] - lower[e];
614      }
615      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {
616        Supply sum = 0;
617        if (n == s) sum =  flow_value;
618        if (n == t) sum = -flow_value;
619        for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)
620          sum += lower[e];
621        for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)
622          sum -= lower[e];
623        _supply[_node_ref[n]] = sum;
624      }
625      _valid_supply = true;
626    }
627
628    /// \brief Simple constructor of the class (without lower bounds).
629    ///
630    /// Simple constructor of the class (without lower bounds).
631    ///
632    /// \param graph The directed graph the algorithm runs on.
633    /// \param capacity The capacities (upper bounds) of the edges.
634    /// \param cost The cost (length) values of the edges.
635    /// \param s The source node.
636    /// \param t The target node.
637    /// \param flow_value The required amount of flow from node \c s
638    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
639    NetworkSimplex( const Graph &graph,
640                    const CapacityMap &capacity,
641                    const CostMap &cost,
642                    typename Graph::Node s,
643                    typename Graph::Node t,
644                    typename SupplyMap::Value flow_value ) :
645      _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
646      _cost(_graph), _supply(_graph, 0), _flow(_graph),
647      _potential(_graph), _depth(_graph), _parent(_graph),
648      _pred_edge(_graph), _thread(_graph), _forward(_graph),
649      _state(_graph), _red_cost(_graph, _cost, _potential),
650      _flow_result(graph), _potential_result(graph),
651      _node_ref(graph), _edge_ref(graph)
652    {
653      // Copying _graph_ref to graph
654      copyGraph(_graph, _graph_ref)
655        .edgeMap(_capacity, capacity)
656        .edgeMap(_cost, cost)
657        .nodeRef(_node_ref)
658        .edgeRef(_edge_ref)
659        .run();
660      _supply[_node_ref[s]] =  flow_value;
661      _supply[_node_ref[t]] = -flow_value;
662      _valid_supply = true;
663    }
664
665    /// \brief Runs the algorithm.
666    ///
667    /// Runs the algorithm.
668    ///
669    /// \param pivot_rule The pivot rule that is used during the
670    /// algorithm.
671    ///
672    /// The available pivot rules:
673    ///
674    /// - FIRST_ELIGIBLE_PIVOT The next eligible edge is selected in
675    /// a wraparound fashion in every iteration
676    /// (\ref FirstEligiblePivotRule).
677    ///
678    /// - BEST_ELIGIBLE_PIVOT The best eligible edge is selected in
679    /// every iteration (\ref BestEligiblePivotRule).
680    ///
681    /// - BLOCK_SEARCH_PIVOT A specified number of edges are examined in
682    /// every iteration in a wraparound fashion and the best eligible
683    /// edge is selected from this block (\ref BlockSearchPivotRule).
684    ///
685    /// - LIMITED_SEARCH_PIVOT A specified number of eligible edges are
686    /// examined in every iteration in a wraparound fashion and the best
687    /// one is selected from them (\ref LimitedSearchPivotRule).
688    ///
689    /// - CANDIDATE_LIST_PIVOT In major iterations a candidate list is
690    /// built from eligible edges and it is used for edge selection in
691    /// the following minor iterations (\ref CandidateListPivotRule).
692    ///
693    /// - COMBINED_PIVOT This is a combined version of the two fastest
694    /// pivot rules.
695    /// For rather sparse graphs \ref LimitedSearchPivotRule
696    /// "Limited Search" implementation is used, otherwise
697    /// \ref BlockSearchPivotRule "Block Search" pivot rule is used.
698    /// According to our benchmark tests this combined method is the
699    /// most efficient.
700    ///
701    /// \return \c true if a feasible flow can be found.
702    bool run(PivotRuleEnum pivot_rule = COMBINED_PIVOT) {
703      return init() && start(pivot_rule);
704    }
705
706    /// \brief Returns a const reference to the edge map storing the
707    /// found flow.
708    ///
709    /// Returns a const reference to the edge map storing the found flow.
710    ///
711    /// \pre \ref run() must be called before using this function.
712    const FlowMap& flowMap() const {
713      return _flow_result;
714    }
715
716    /// \brief Returns a const reference to the node map storing the
717    /// found potentials (the dual solution).
718    ///
719    /// Returns a const reference to the node map storing the found
720    /// potentials (the dual solution).
721    ///
722    /// \pre \ref run() must be called before using this function.
723    const PotentialMap& potentialMap() const {
724      return _potential_result;
725    }
726
727    /// \brief Returns the total cost of the found flow.
728    ///
729    /// Returns the total cost of the found flow. The complexity of the
730    /// function is \f$ O(e) \f$.
731    ///
732    /// \pre \ref run() must be called before using this function.
733    Cost totalCost() const {
734      Cost c = 0;
735      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
736        c += _flow_result[e] * _cost[_edge_ref[e]];
737      return c;
738    }
739
740  private:
741
742    /// \brief Extends the underlaying graph and initializes all the
743    /// node and edge maps.
744    bool init() {
745      if (!_valid_supply) return false;
746
747      // Initializing state and flow maps
748      for (EdgeIt e(_graph); e != INVALID; ++e) {
749        _flow[e] = 0;
750        _state[e] = STATE_LOWER;
751      }
752
753      // Adding an artificial root node to the graph
754      _root = _graph.addNode();
755      _parent[_root] = INVALID;
756      _pred_edge[_root] = INVALID;
757      _depth[_root] = 0;
758      _supply[_root] = 0;
759      _potential[_root] = 0;
760
761      // Adding artificial edges to the graph and initializing the node
762      // maps of the spanning tree data structure
763      Node last = _root;
764      Edge e;
765      Cost max_cost = std::numeric_limits<Cost>::max() / 4;
766      for (NodeIt u(_graph); u != INVALID; ++u) {
767        if (u == _root) continue;
768        _thread[last] = u;
769        last = u;
770        _parent[u] = _root;
771        _depth[u] = 1;
772        if (_supply[u] >= 0) {
773          e = _graph.addEdge(u, _root);
774          _flow[e] = _supply[u];
775          _forward[u] = true;
776          _potential[u] = -max_cost;
777        } else {
778          e = _graph.addEdge(_root, u);
779          _flow[e] = -_supply[u];
780          _forward[u] = false;
781          _potential[u] = max_cost;
782        }
783        _cost[e] = max_cost;
784        _capacity[e] = std::numeric_limits<Capacity>::max();
785        _state[e] = STATE_TREE;
786        _pred_edge[u] = e;
787      }
788      _thread[last] = _root;
789
790      return true;
791    }
792
793    /// Finds the join node.
794    Node findJoinNode() {
795      Node u = _graph.source(_in_edge);
796      Node v = _graph.target(_in_edge);
797      while (u != v) {
798        if (_depth[u] == _depth[v]) {
799          u = _parent[u];
800          v = _parent[v];
801        }
802        else if (_depth[u] > _depth[v]) u = _parent[u];
803        else v = _parent[v];
804      }
805      return u;
806    }
807
808    /// \brief Finds the leaving edge of the cycle. Returns \c true if
809    /// the leaving edge is not the same as the entering edge.
810    bool findLeavingEdge() {
811      // Initializing first and second nodes according to the direction
812      // of the cycle
813      if (_state[_in_edge] == STATE_LOWER) {
814        first  = _graph.source(_in_edge);
815        second = _graph.target(_in_edge);
816      } else {
817        first  = _graph.target(_in_edge);
818        second = _graph.source(_in_edge);
819      }
820      delta = _capacity[_in_edge];
821      bool result = false;
822      Capacity d;
823      Edge e;
824
825      // Searching the cycle along the path form the first node to the
826      // root node
827      for (Node u = first; u != join; u = _parent[u]) {
828        e = _pred_edge[u];
829        d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e];
830        if (d < delta) {
831          delta = d;
832          u_out = u;
833          u_in = first;
834          v_in = second;
835          result = true;
836        }
837      }
838      // Searching the cycle along the path form the second node to the
839      // root node
840      for (Node u = second; u != join; u = _parent[u]) {
841        e = _pred_edge[u];
842        d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e];
843        if (d <= delta) {
844          delta = d;
845          u_out = u;
846          u_in = second;
847          v_in = first;
848          result = true;
849        }
850      }
851      return result;
852    }
853
854    /// Changes \c flow and \c state edge maps.
855    void changeFlows(bool change) {
856      // Augmenting along the cycle
857      if (delta > 0) {
858        Capacity val = _state[_in_edge] * delta;
859        _flow[_in_edge] += val;
860        for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) {
861          _flow[_pred_edge[u]] += _forward[u] ? -val : val;
862        }
863        for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) {
864          _flow[_pred_edge[u]] += _forward[u] ? val : -val;
865        }
866      }
867      // Updating the state of the entering and leaving edges
868      if (change) {
869        _state[_in_edge] = STATE_TREE;
870        _state[_pred_edge[u_out]] =
871          (_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
872      } else {
873        _state[_in_edge] = -_state[_in_edge];
874      }
875    }
876
877    /// Updates \c thread and \c parent node maps.
878    void updateThreadParent() {
879      Node u;
880      v_out = _parent[u_out];
881
882      // Handling the case when join and v_out coincide
883      bool par_first = false;
884      if (join == v_out) {
885        for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
886        if (u == v_in) {
887          par_first = true;
888          while (_thread[u] != u_out) u = _thread[u];
889          first = u;
890        }
891      }
892
893      // Finding the last successor of u_in (u) and the node after it
894      // (right) according to the thread index
895      for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
896      right = _thread[u];
897      if (_thread[v_in] == u_out) {
898        for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;
899        if (last == u_out) last = _thread[last];
900      }
901      else last = _thread[v_in];
902
903      // Updating stem nodes
904      _thread[v_in] = stem = u_in;
905      par_stem = v_in;
906      while (stem != u_out) {
907        _thread[u] = new_stem = _parent[stem];
908
909        // Finding the node just before the stem node (u) according to
910        // the original thread index
911        for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
912        _thread[u] = right;
913
914        // Changing the parent node of stem and shifting stem and
915        // par_stem nodes
916        _parent[stem] = par_stem;
917        par_stem = stem;
918        stem = new_stem;
919
920        // Finding the last successor of stem (u) and the node after it
921        // (right) according to the thread index
922        for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
923        right = _thread[u];
924      }
925      _parent[u_out] = par_stem;
926      _thread[u] = last;
927
928      if (join == v_out && par_first) {
929        if (first != v_in) _thread[first] = right;
930      } else {
931        for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
932        _thread[u] = right;
933      }
934    }
935
936    /// Updates \c pred_edge and \c forward node maps.
937    void updatePredEdge() {
938      Node u = u_out, v;
939      while (u != u_in) {
940        v = _parent[u];
941        _pred_edge[u] = _pred_edge[v];
942        _forward[u] = !_forward[v];
943        u = v;
944      }
945      _pred_edge[u_in] = _in_edge;
946      _forward[u_in] = (u_in == _graph.source(_in_edge));
947    }
948
949    /// Updates \c depth and \c potential node maps.
950    void updateDepthPotential() {
951      _depth[u_in] = _depth[v_in] + 1;
952      _potential[u_in] = _forward[u_in] ?
953        _potential[v_in] - _cost[_pred_edge[u_in]] :
954        _potential[v_in] + _cost[_pred_edge[u_in]];
955
956      Node u = _thread[u_in], v;
957      while (true) {
958        v = _parent[u];
959        if (v == INVALID) break;
960        _depth[u] = _depth[v] + 1;
961        _potential[u] = _forward[u] ?
962          _potential[v] - _cost[_pred_edge[u]] :
963          _potential[v] + _cost[_pred_edge[u]];
964        if (_depth[u] <= _depth[v_in]) break;
965        u = _thread[u];
966      }
967    }
968
969    /// Executes the algorithm.
970    bool start(PivotRuleEnum pivot_rule) {
971      switch (pivot_rule) {
972        case FIRST_ELIGIBLE_PIVOT:
973          return start<FirstEligiblePivotRule>();
974        case BEST_ELIGIBLE_PIVOT:
975          return start<BestEligiblePivotRule>();
976        case BLOCK_SEARCH_PIVOT:
977          return start<BlockSearchPivotRule>();
978        case LIMITED_SEARCH_PIVOT:
979          return start<LimitedSearchPivotRule>();
980        case CANDIDATE_LIST_PIVOT:
981          return start<CandidateListPivotRule>();
982        case COMBINED_PIVOT:
983          if ( countEdges(_graph) / countNodes(_graph) <=
984               COMBINED_PIVOT_MAX_DEG )
985            return start<LimitedSearchPivotRule>();
986          else
987            return start<BlockSearchPivotRule>();
988      }
989      return false;
990    }
991
992    template<class PivotRuleImplementation>
993    bool start() {
994      PivotRuleImplementation pivot(*this);
995
996      // Executing the network simplex algorithm
997      while (pivot.findEnteringEdge()) {
998        join = findJoinNode();
999        bool change = findLeavingEdge();
1000        changeFlows(change);
1001        if (change) {
1002          updateThreadParent();
1003          updatePredEdge();
1004          updateDepthPotential();
1005        }
1006      }
1007
1008      // Checking if the flow amount equals zero on all the artificial
1009      // edges
1010      for (InEdgeIt e(_graph, _root); e != INVALID; ++e)
1011        if (_flow[e] > 0) return false;
1012      for (OutEdgeIt e(_graph, _root); e != INVALID; ++e)
1013        if (_flow[e] > 0) return false;
1014
1015      // Copying flow values to _flow_result
1016      if (_lower) {
1017        for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
1018          _flow_result[e] = (*_lower)[e] + _flow[_edge_ref[e]];
1019      } else {
1020        for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
1021          _flow_result[e] = _flow[_edge_ref[e]];
1022      }
1023      // Copying potential values to _potential_result
1024      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
1025        _potential_result[n] = _potential[_node_ref[n]];
1026
1027      return true;
1028    }
1029
1030  }; //class NetworkSimplex
1031
1032  ///@}
1033
1034} //namespace lemon
1035
1036#endif //LEMON_NETWORK_SIMPLEX_H
Note: See TracBrowser for help on using the repository browser.