COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/network_simplex.h @ 2629:84354c78b068

Last change on this file since 2629:84354c78b068 was 2629:84354c78b068, checked in by Peter Kovacs, 12 years ago

Improved constructors for min cost flow classes
Removing the non-zero lower bounds is faster

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