COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/network_simplex.h @ 2587:061be2e64eca

Last change on this file since 2587:061be2e64eca was 2581:054566ac0934, checked in by Peter Kovacs, 16 years ago

Query improvements in the min cost flow algorithms.

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