COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/network_simplex.h @ 2593:8eed667ea23c

Last change on this file since 2593:8eed667ea23c was 2593:8eed667ea23c, checked in by Peter Kovacs, 17 years ago

Fix static member initializations (ticket #30).

File size: 35.3 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2008
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
12 *
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_NETWORK_SIMPLEX_H
20#define LEMON_NETWORK_SIMPLEX_H
21
22/// \ingroup min_cost_flow
23///
24/// \file
25/// \brief Network simplex algorithm for finding a minimum cost flow.
26
27#include <vector>
28#include <limits>
29
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 SAMPLE_SIZE_FACTOR = 15;
284      static const int MIN_SAMPLE_SIZE = 10;
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 = countEdges(_ns._graph) *
293                       SAMPLE_SIZE_FACTOR / 10000;
294        if (_sample_size < MIN_SAMPLE_SIZE)
295          _sample_size = MIN_SAMPLE_SIZE;
296      }
297
298      /// Finds the next entering edge.
299      bool findEnteringEdge() {
300        Cost curr, min = 0;
301        int cnt = 0;
302        for (EdgeIt e = _next_edge; e != INVALID; ++e) {
303          if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
304            min = curr;
305            _min_edge = e;
306          }
307          if (curr < 0 && ++cnt == _sample_size) break;
308        }
309        if (min == 0) {
310          for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
311            if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
312              min = curr;
313              _min_edge = e;
314            }
315            if (curr < 0 && ++cnt == _sample_size) break;
316          }
317        }
318        _ns._in_edge = _min_edge;
319        _next_edge = ++_min_edge;
320        return min < 0;
321      }
322    }; //class LimitedSearchPivotRule
323
324    /// \brief Implementation of the "Candidate List" pivot rule for the
325    /// \ref NetworkSimplex "network simplex" algorithm.
326    ///
327    /// This class implements the "Candidate List" pivot rule
328    /// for the \ref NetworkSimplex "network simplex" algorithm.
329    class CandidateListPivotRule
330    {
331    private:
332
333      NetworkSimplex &_ns;
334
335      // The list of candidate edges.
336      std::vector<Edge> _candidates;
337      // The maximum length of the edge list.
338      int _list_length;
339      // The maximum number of minor iterations between two major
340      // itarations.
341      int _minor_limit;
342
343      int _minor_count;
344      EdgeIt _next_edge;
345
346      static const int LIST_LENGTH_FACTOR = 20;
347      static const int MINOR_LIMIT_FACTOR = 10;
348      static const int MIN_LIST_LENGTH = 10;
349      static const int MIN_MINOR_LIMIT = 2;
350
351    public:
352
353      /// Constructor.
354      CandidateListPivotRule(NetworkSimplex &ns) :
355        _ns(ns), _next_edge(ns._graph)
356      {
357        int edge_num = countEdges(_ns._graph);
358        _minor_count = 0;
359        _list_length = edge_num * LIST_LENGTH_FACTOR / 10000;
360        if (_list_length < MIN_LIST_LENGTH)
361          _list_length = MIN_LIST_LENGTH;
362        _minor_limit = _list_length * MINOR_LIMIT_FACTOR / 100;
363        if (_minor_limit < MIN_MINOR_LIMIT)
364          _minor_limit = MIN_MINOR_LIMIT;
365      }
366
367      /// Finds the next entering edge.
368      bool findEnteringEdge() {
369        Cost min, curr;
370        if (_minor_count < _minor_limit && _candidates.size() > 0) {
371          // Minor iteration
372          ++_minor_count;
373          Edge e;
374          min = 0;
375          for (int i = 0; i < int(_candidates.size()); ++i) {
376            e = _candidates[i];
377            if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
378              min = curr;
379              _ns._in_edge = e;
380            }
381          }
382          if (min < 0) return true;
383        }
384
385        // Major iteration
386        _candidates.clear();
387        EdgeIt e = _next_edge;
388        min = 0;
389        for ( ; e != INVALID; ++e) {
390          if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
391            _candidates.push_back(e);
392            if (curr < min) {
393              min = curr;
394              _ns._in_edge = e;
395            }
396            if (int(_candidates.size()) == _list_length) break;
397          }
398        }
399        if (int(_candidates.size()) < _list_length) {
400          for (e = EdgeIt(_ns._graph); e != _next_edge; ++e) {
401            if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
402              _candidates.push_back(e);
403              if (curr < min) {
404                min = curr;
405                _ns._in_edge = e;
406              }
407              if (int(_candidates.size()) == _list_length) break;
408            }
409          }
410        }
411        if (_candidates.size() == 0) return false;
412        _minor_count = 1;
413        _next_edge = ++e;
414        return true;
415      }
416    }; //class CandidateListPivotRule
417
418  private:
419
420    // State constants for edges
421    enum EdgeStateEnum {
422      STATE_UPPER = -1,
423      STATE_TREE  =  0,
424      STATE_LOWER =  1
425    };
426
427    // Constant for the combined pivot rule.
428    static const int COMBINED_PIVOT_MAX_DEG = 5;
429
430  private:
431
432    // The directed graph the algorithm runs on
433    SGraph _graph;
434    // The original graph
435    const Graph &_graph_ref;
436    // The original lower bound map
437    const LowerMap *_lower;
438    // The capacity map
439    SCapacityMap _capacity;
440    // The cost map
441    SCostMap _cost;
442    // The supply map
443    SSupplyMap _supply;
444    bool _valid_supply;
445
446    // Edge map of the current flow
447    SCapacityMap _flow;
448    // Node map of the current potentials
449    SPotentialMap _potential;
450
451    // The depth node map of the spanning tree structure
452    IntNodeMap _depth;
453    // The parent node map of the spanning tree structure
454    NodeNodeMap _parent;
455    // The pred_edge node map of the spanning tree structure
456    EdgeNodeMap _pred_edge;
457    // The thread node map of the spanning tree structure
458    NodeNodeMap _thread;
459    // The forward node map of the spanning tree structure
460    BoolNodeMap _forward;
461    // The state edge map
462    IntEdgeMap _state;
463    // The root node of the starting spanning tree
464    Node _root;
465
466    // The reduced cost map
467    ReducedCostMap _red_cost;
468
469    // Members for handling the original graph
470    FlowMap *_flow_result;
471    PotentialMap *_potential_result;
472    bool _local_flow;
473    bool _local_potential;
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 (with lower bounds).
491    ///
492    /// General constructor (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(0), _potential_result(0),
510      _local_flow(false), _local_potential(false),
511      _node_ref(graph), _edge_ref(graph)
512    {
513      // Checking the sum of supply values
514      Supply sum = 0;
515      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
516        sum += supply[n];
517      if (!(_valid_supply = sum == 0)) return;
518
519      // Copying _graph_ref to _graph
520      _graph.reserveNode(countNodes(_graph_ref) + 1);
521      _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
522      copyGraph(_graph, _graph_ref)
523        .edgeMap(_cost, cost)
524        .nodeRef(_node_ref)
525        .edgeRef(_edge_ref)
526        .run();
527
528      // Removing non-zero lower bounds
529      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
530        _capacity[_edge_ref[e]] = capacity[e] - lower[e];
531      }
532      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {
533        Supply s = supply[n];
534        for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)
535          s += lower[e];
536        for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)
537          s -= lower[e];
538        _supply[_node_ref[n]] = s;
539      }
540    }
541
542    /// \brief General constructor (without lower bounds).
543    ///
544    /// General constructor (without lower bounds).
545    ///
546    /// \param graph The directed graph the algorithm runs on.
547    /// \param capacity The capacities (upper bounds) of the edges.
548    /// \param cost The cost (length) values of the edges.
549    /// \param supply The supply values of the nodes (signed).
550    NetworkSimplex( const Graph &graph,
551                    const CapacityMap &capacity,
552                    const CostMap &cost,
553                    const SupplyMap &supply ) :
554      _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
555      _cost(_graph), _supply(_graph), _flow(_graph),
556      _potential(_graph), _depth(_graph), _parent(_graph),
557      _pred_edge(_graph), _thread(_graph), _forward(_graph),
558      _state(_graph), _red_cost(_graph, _cost, _potential),
559      _flow_result(0), _potential_result(0),
560      _local_flow(false), _local_potential(false),
561      _node_ref(graph), _edge_ref(graph)
562    {
563      // Checking the sum of supply values
564      Supply sum = 0;
565      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
566        sum += supply[n];
567      if (!(_valid_supply = sum == 0)) return;
568
569      // Copying _graph_ref to graph
570      copyGraph(_graph, _graph_ref)
571        .edgeMap(_capacity, capacity)
572        .edgeMap(_cost, cost)
573        .nodeMap(_supply, supply)
574        .nodeRef(_node_ref)
575        .edgeRef(_edge_ref)
576        .run();
577    }
578
579    /// \brief Simple constructor (with lower bounds).
580    ///
581    /// Simple constructor (with lower bounds).
582    ///
583    /// \param graph The directed graph the algorithm runs on.
584    /// \param lower The lower bounds of the edges.
585    /// \param capacity The capacities (upper bounds) of the edges.
586    /// \param cost The cost (length) values of the edges.
587    /// \param s The source node.
588    /// \param t The target node.
589    /// \param flow_value The required amount of flow from node \c s
590    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
591    NetworkSimplex( const Graph &graph,
592                    const LowerMap &lower,
593                    const CapacityMap &capacity,
594                    const CostMap &cost,
595                    typename Graph::Node s,
596                    typename Graph::Node t,
597                    typename SupplyMap::Value flow_value ) :
598      _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
599      _cost(_graph), _supply(_graph), _flow(_graph),
600      _potential(_graph), _depth(_graph), _parent(_graph),
601      _pred_edge(_graph), _thread(_graph), _forward(_graph),
602      _state(_graph), _red_cost(_graph, _cost, _potential),
603      _flow_result(0), _potential_result(0),
604      _local_flow(false), _local_potential(false),
605      _node_ref(graph), _edge_ref(graph)
606    {
607      // Copying _graph_ref to graph
608      copyGraph(_graph, _graph_ref)
609        .edgeMap(_cost, cost)
610        .nodeRef(_node_ref)
611        .edgeRef(_edge_ref)
612        .run();
613
614      // Removing non-zero lower bounds
615      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
616        _capacity[_edge_ref[e]] = capacity[e] - lower[e];
617      }
618      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {
619        Supply sum = 0;
620        if (n == s) sum =  flow_value;
621        if (n == t) sum = -flow_value;
622        for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)
623          sum += lower[e];
624        for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)
625          sum -= lower[e];
626        _supply[_node_ref[n]] = sum;
627      }
628      _valid_supply = true;
629    }
630
631    /// \brief Simple constructor (without lower bounds).
632    ///
633    /// Simple constructor (without lower bounds).
634    ///
635    /// \param graph The directed graph the algorithm runs on.
636    /// \param capacity The capacities (upper bounds) of the edges.
637    /// \param cost The cost (length) values of the edges.
638    /// \param s The source node.
639    /// \param t The target node.
640    /// \param flow_value The required amount of flow from node \c s
641    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
642    NetworkSimplex( const Graph &graph,
643                    const CapacityMap &capacity,
644                    const CostMap &cost,
645                    typename Graph::Node s,
646                    typename Graph::Node t,
647                    typename SupplyMap::Value flow_value ) :
648      _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
649      _cost(_graph), _supply(_graph, 0), _flow(_graph),
650      _potential(_graph), _depth(_graph), _parent(_graph),
651      _pred_edge(_graph), _thread(_graph), _forward(_graph),
652      _state(_graph), _red_cost(_graph, _cost, _potential),
653      _flow_result(0), _potential_result(0),
654      _local_flow(false), _local_potential(false),
655      _node_ref(graph), _edge_ref(graph)
656    {
657      // Copying _graph_ref to graph
658      copyGraph(_graph, _graph_ref)
659        .edgeMap(_capacity, capacity)
660        .edgeMap(_cost, cost)
661        .nodeRef(_node_ref)
662        .edgeRef(_edge_ref)
663        .run();
664      _supply[_node_ref[s]] =  flow_value;
665      _supply[_node_ref[t]] = -flow_value;
666      _valid_supply = true;
667    }
668
669    /// Destructor.
670    ~NetworkSimplex() {
671      if (_local_flow) delete _flow_result;
672      if (_local_potential) delete _potential_result;
673    }
674
675    /// \brief Sets the flow map.
676    ///
677    /// Sets the flow map.
678    ///
679    /// \return \c (*this)
680    NetworkSimplex& flowMap(FlowMap &map) {
681      if (_local_flow) {
682        delete _flow_result;
683        _local_flow = false;
684      }
685      _flow_result = &map;
686      return *this;
687    }
688
689    /// \brief Sets the potential map.
690    ///
691    /// Sets the potential map.
692    ///
693    /// \return \c (*this)
694    NetworkSimplex& potentialMap(PotentialMap &map) {
695      if (_local_potential) {
696        delete _potential_result;
697        _local_potential = false;
698      }
699      _potential_result = &map;
700      return *this;
701    }
702
703    /// \name Execution control
704    /// The only way to execute the algorithm is to call the run()
705    /// function.
706
707    /// @{
708
709    /// \brief Runs the algorithm.
710    ///
711    /// Runs the algorithm.
712    ///
713    /// \param pivot_rule The pivot rule that is used during the
714    /// algorithm.
715    ///
716    /// The available pivot rules:
717    ///
718    /// - FIRST_ELIGIBLE_PIVOT The next eligible edge is selected in
719    /// a wraparound fashion in every iteration
720    /// (\ref FirstEligiblePivotRule).
721    ///
722    /// - BEST_ELIGIBLE_PIVOT The best eligible edge is selected in
723    /// every iteration (\ref BestEligiblePivotRule).
724    ///
725    /// - BLOCK_SEARCH_PIVOT A specified number of edges are examined in
726    /// every iteration in a wraparound fashion and the best eligible
727    /// edge is selected from this block (\ref BlockSearchPivotRule).
728    ///
729    /// - LIMITED_SEARCH_PIVOT A specified number of eligible edges are
730    /// examined in every iteration in a wraparound fashion and the best
731    /// one is selected from them (\ref LimitedSearchPivotRule).
732    ///
733    /// - CANDIDATE_LIST_PIVOT In major iterations a candidate list is
734    /// built from eligible edges and it is used for edge selection in
735    /// the following minor iterations (\ref CandidateListPivotRule).
736    ///
737    /// - COMBINED_PIVOT This is a combined version of the two fastest
738    /// pivot rules.
739    /// For rather sparse graphs \ref LimitedSearchPivotRule
740    /// "Limited Search" implementation is used, otherwise
741    /// \ref BlockSearchPivotRule "Block Search" pivot rule is used.
742    /// According to our benchmark tests this combined method is the
743    /// most efficient.
744    ///
745    /// \return \c true if a feasible flow can be found.
746    bool run(PivotRuleEnum pivot_rule = COMBINED_PIVOT) {
747      return init() && start(pivot_rule);
748    }
749
750    /// @}
751
752    /// \name Query Functions
753    /// The result of the algorithm can be obtained using these
754    /// functions.
755    /// \n run() must be called before using them.
756
757    /// @{
758
759    /// \brief Returns a const reference to the edge map storing the
760    /// found flow.
761    ///
762    /// Returns a const reference to the edge map storing the found flow.
763    ///
764    /// \pre \ref run() must be called before using this function.
765    const FlowMap& flowMap() const {
766      return *_flow_result;
767    }
768
769    /// \brief Returns a const reference to the node map storing the
770    /// found potentials (the dual solution).
771    ///
772    /// Returns a const reference to the node map storing the found
773    /// potentials (the dual solution).
774    ///
775    /// \pre \ref run() must be called before using this function.
776    const PotentialMap& potentialMap() const {
777      return *_potential_result;
778    }
779
780    /// \brief Returns the flow on the given edge.
781    ///
782    /// Returns the flow on the given edge.
783    ///
784    /// \pre \ref run() must be called before using this function.
785    Capacity flow(const typename Graph::Edge& edge) const {
786      return (*_flow_result)[edge];
787    }
788
789    /// \brief Returns the potential of the given node.
790    ///
791    /// Returns the potential of the given node.
792    ///
793    /// \pre \ref run() must be called before using this function.
794    Cost potential(const typename Graph::Node& node) const {
795      return (*_potential_result)[node];
796    }
797
798    /// \brief Returns the total cost of the found flow.
799    ///
800    /// Returns the total cost of the found flow. The complexity of the
801    /// function is \f$ O(e) \f$.
802    ///
803    /// \pre \ref run() must be called before using this function.
804    Cost totalCost() const {
805      Cost c = 0;
806      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
807        c += (*_flow_result)[e] * _cost[_edge_ref[e]];
808      return c;
809    }
810
811    /// @}
812
813  private:
814
815    /// \brief Extends the underlaying graph and initializes all the
816    /// node and edge maps.
817    bool init() {
818      if (!_valid_supply) return false;
819
820      // Initializing result maps
821      if (!_flow_result) {
822        _flow_result = new FlowMap(_graph_ref);
823        _local_flow = true;
824      }
825      if (!_potential_result) {
826        _potential_result = new PotentialMap(_graph_ref);
827        _local_potential = true;
828      }
829
830      // Initializing state and flow maps
831      for (EdgeIt e(_graph); e != INVALID; ++e) {
832        _flow[e] = 0;
833        _state[e] = STATE_LOWER;
834      }
835
836      // Adding an artificial root node to the graph
837      _root = _graph.addNode();
838      _parent[_root] = INVALID;
839      _pred_edge[_root] = INVALID;
840      _depth[_root] = 0;
841      _supply[_root] = 0;
842      _potential[_root] = 0;
843
844      // Adding artificial edges to the graph and initializing the node
845      // maps of the spanning tree data structure
846      Node last = _root;
847      Edge e;
848      Cost max_cost = std::numeric_limits<Cost>::max() / 4;
849      for (NodeIt u(_graph); u != INVALID; ++u) {
850        if (u == _root) continue;
851        _thread[last] = u;
852        last = u;
853        _parent[u] = _root;
854        _depth[u] = 1;
855        if (_supply[u] >= 0) {
856          e = _graph.addEdge(u, _root);
857          _flow[e] = _supply[u];
858          _forward[u] = true;
859          _potential[u] = -max_cost;
860        } else {
861          e = _graph.addEdge(_root, u);
862          _flow[e] = -_supply[u];
863          _forward[u] = false;
864          _potential[u] = max_cost;
865        }
866        _cost[e] = max_cost;
867        _capacity[e] = std::numeric_limits<Capacity>::max();
868        _state[e] = STATE_TREE;
869        _pred_edge[u] = e;
870      }
871      _thread[last] = _root;
872
873      return true;
874    }
875
876    /// Finds the join node.
877    Node findJoinNode() {
878      Node u = _graph.source(_in_edge);
879      Node v = _graph.target(_in_edge);
880      while (u != v) {
881        if (_depth[u] == _depth[v]) {
882          u = _parent[u];
883          v = _parent[v];
884        }
885        else if (_depth[u] > _depth[v]) u = _parent[u];
886        else v = _parent[v];
887      }
888      return u;
889    }
890
891    /// \brief Finds the leaving edge of the cycle. Returns \c true if
892    /// the leaving edge is not the same as the entering edge.
893    bool findLeavingEdge() {
894      // Initializing first and second nodes according to the direction
895      // of the cycle
896      if (_state[_in_edge] == STATE_LOWER) {
897        first  = _graph.source(_in_edge);
898        second = _graph.target(_in_edge);
899      } else {
900        first  = _graph.target(_in_edge);
901        second = _graph.source(_in_edge);
902      }
903      delta = _capacity[_in_edge];
904      bool result = false;
905      Capacity d;
906      Edge e;
907
908      // Searching the cycle along the path form the first node to the
909      // root node
910      for (Node u = first; u != join; u = _parent[u]) {
911        e = _pred_edge[u];
912        d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e];
913        if (d < delta) {
914          delta = d;
915          u_out = u;
916          u_in = first;
917          v_in = second;
918          result = true;
919        }
920      }
921      // Searching the cycle along the path form the second node to the
922      // root node
923      for (Node u = second; u != join; u = _parent[u]) {
924        e = _pred_edge[u];
925        d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e];
926        if (d <= delta) {
927          delta = d;
928          u_out = u;
929          u_in = second;
930          v_in = first;
931          result = true;
932        }
933      }
934      return result;
935    }
936
937    /// Changes \c flow and \c state edge maps.
938    void changeFlows(bool change) {
939      // Augmenting along the cycle
940      if (delta > 0) {
941        Capacity val = _state[_in_edge] * delta;
942        _flow[_in_edge] += val;
943        for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) {
944          _flow[_pred_edge[u]] += _forward[u] ? -val : val;
945        }
946        for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) {
947          _flow[_pred_edge[u]] += _forward[u] ? val : -val;
948        }
949      }
950      // Updating the state of the entering and leaving edges
951      if (change) {
952        _state[_in_edge] = STATE_TREE;
953        _state[_pred_edge[u_out]] =
954          (_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
955      } else {
956        _state[_in_edge] = -_state[_in_edge];
957      }
958    }
959
960    /// Updates \c thread and \c parent node maps.
961    void updateThreadParent() {
962      Node u;
963      v_out = _parent[u_out];
964
965      // Handling the case when join and v_out coincide
966      bool par_first = false;
967      if (join == v_out) {
968        for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
969        if (u == v_in) {
970          par_first = true;
971          while (_thread[u] != u_out) u = _thread[u];
972          first = u;
973        }
974      }
975
976      // Finding the last successor of u_in (u) and the node after it
977      // (right) according to the thread index
978      for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
979      right = _thread[u];
980      if (_thread[v_in] == u_out) {
981        for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;
982        if (last == u_out) last = _thread[last];
983      }
984      else last = _thread[v_in];
985
986      // Updating stem nodes
987      _thread[v_in] = stem = u_in;
988      par_stem = v_in;
989      while (stem != u_out) {
990        _thread[u] = new_stem = _parent[stem];
991
992        // Finding the node just before the stem node (u) according to
993        // the original thread index
994        for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
995        _thread[u] = right;
996
997        // Changing the parent node of stem and shifting stem and
998        // par_stem nodes
999        _parent[stem] = par_stem;
1000        par_stem = stem;
1001        stem = new_stem;
1002
1003        // Finding the last successor of stem (u) and the node after it
1004        // (right) according to the thread index
1005        for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
1006        right = _thread[u];
1007      }
1008      _parent[u_out] = par_stem;
1009      _thread[u] = last;
1010
1011      if (join == v_out && par_first) {
1012        if (first != v_in) _thread[first] = right;
1013      } else {
1014        for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
1015        _thread[u] = right;
1016      }
1017    }
1018
1019    /// Updates \c pred_edge and \c forward node maps.
1020    void updatePredEdge() {
1021      Node u = u_out, v;
1022      while (u != u_in) {
1023        v = _parent[u];
1024        _pred_edge[u] = _pred_edge[v];
1025        _forward[u] = !_forward[v];
1026        u = v;
1027      }
1028      _pred_edge[u_in] = _in_edge;
1029      _forward[u_in] = (u_in == _graph.source(_in_edge));
1030    }
1031
1032    /// Updates \c depth and \c potential node maps.
1033    void updateDepthPotential() {
1034      _depth[u_in] = _depth[v_in] + 1;
1035      _potential[u_in] = _forward[u_in] ?
1036        _potential[v_in] - _cost[_pred_edge[u_in]] :
1037        _potential[v_in] + _cost[_pred_edge[u_in]];
1038
1039      Node u = _thread[u_in], v;
1040      while (true) {
1041        v = _parent[u];
1042        if (v == INVALID) break;
1043        _depth[u] = _depth[v] + 1;
1044        _potential[u] = _forward[u] ?
1045          _potential[v] - _cost[_pred_edge[u]] :
1046          _potential[v] + _cost[_pred_edge[u]];
1047        if (_depth[u] <= _depth[v_in]) break;
1048        u = _thread[u];
1049      }
1050    }
1051
1052    /// Executes the algorithm.
1053    bool start(PivotRuleEnum pivot_rule) {
1054      switch (pivot_rule) {
1055        case FIRST_ELIGIBLE_PIVOT:
1056          return start<FirstEligiblePivotRule>();
1057        case BEST_ELIGIBLE_PIVOT:
1058          return start<BestEligiblePivotRule>();
1059        case BLOCK_SEARCH_PIVOT:
1060          return start<BlockSearchPivotRule>();
1061        case LIMITED_SEARCH_PIVOT:
1062          return start<LimitedSearchPivotRule>();
1063        case CANDIDATE_LIST_PIVOT:
1064          return start<CandidateListPivotRule>();
1065        case COMBINED_PIVOT:
1066          if ( countEdges(_graph) / countNodes(_graph) <=
1067               COMBINED_PIVOT_MAX_DEG )
1068            return start<LimitedSearchPivotRule>();
1069          else
1070            return start<BlockSearchPivotRule>();
1071      }
1072      return false;
1073    }
1074
1075    template<class PivotRuleImplementation>
1076    bool start() {
1077      PivotRuleImplementation pivot(*this);
1078
1079      // Executing the network simplex algorithm
1080      while (pivot.findEnteringEdge()) {
1081        join = findJoinNode();
1082        bool change = findLeavingEdge();
1083        changeFlows(change);
1084        if (change) {
1085          updateThreadParent();
1086          updatePredEdge();
1087          updateDepthPotential();
1088        }
1089      }
1090
1091      // Checking if the flow amount equals zero on all the artificial
1092      // edges
1093      for (InEdgeIt e(_graph, _root); e != INVALID; ++e)
1094        if (_flow[e] > 0) return false;
1095      for (OutEdgeIt e(_graph, _root); e != INVALID; ++e)
1096        if (_flow[e] > 0) return false;
1097
1098      // Copying flow values to _flow_result
1099      if (_lower) {
1100        for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
1101          (*_flow_result)[e] = (*_lower)[e] + _flow[_edge_ref[e]];
1102      } else {
1103        for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
1104          (*_flow_result)[e] = _flow[_edge_ref[e]];
1105      }
1106      // Copying potential values to _potential_result
1107      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
1108        (*_potential_result)[n] = _potential[_node_ref[n]];
1109
1110      return true;
1111    }
1112
1113  }; //class NetworkSimplex
1114
1115  ///@}
1116
1117} //namespace lemon
1118
1119#endif //LEMON_NETWORK_SIMPLEX_H
Note: See TracBrowser for help on using the repository browser.