COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/network_simplex.h @ 2623:90defb96ee61

Last change on this file since 2623:90defb96ee61 was 2623:90defb96ee61, checked in by Peter Kovacs, 11 years ago

Add missing pointer initializing in min cost flow classes

File size: 38.5 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2008
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
12 *
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_NETWORK_SIMPLEX_H
20#define LEMON_NETWORK_SIMPLEX_H
21
22/// \ingroup min_cost_flow
23///
24/// \file
25/// \brief Network simplex algorithm for finding a minimum cost flow.
26
27#include <vector>
28#include <limits>
29
30#include <lemon/graph_adaptor.h>
31#include <lemon/graph_utils.h>
32#include <lemon/smart_graph.h>
33#include <lemon/math.h>
34
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      // Checking 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      // Copying _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(_cost, cost)
625        .nodeRef(_node_ref)
626        .edgeRef(_edge_ref)
627        .run();
628
629      // Removing non-zero lower bounds
630      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
631        _capacity[_edge_ref[e]] = capacity[e] - lower[e];
632      }
633      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {
634        Supply s = supply[n];
635        for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)
636          s += lower[e];
637        for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)
638          s -= lower[e];
639        _supply[_node_ref[n]] = s;
640      }
641    }
642
643    /// \brief General constructor (without lower bounds).
644    ///
645    /// General constructor (without lower bounds).
646    ///
647    /// \param graph The directed graph the algorithm runs on.
648    /// \param capacity The capacities (upper bounds) of the edges.
649    /// \param cost The cost (length) values of the edges.
650    /// \param supply The supply values of the nodes (signed).
651    NetworkSimplex( const Graph &graph,
652                    const CapacityMap &capacity,
653                    const CostMap &cost,
654                    const SupplyMap &supply ) :
655      _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
656      _cost(_graph), _supply(_graph), _flow(_graph),
657      _potential(_graph), _depth(_graph), _parent(_graph),
658      _pred_edge(_graph), _thread(_graph), _forward(_graph),
659      _state(_graph), _red_cost(_graph, _cost, _potential),
660      _flow_result(NULL), _potential_result(NULL),
661      _local_flow(false), _local_potential(false),
662      _node_ref(graph), _edge_ref(graph)
663    {
664      // Checking the sum of supply values
665      Supply sum = 0;
666      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
667        sum += supply[n];
668      if (!(_valid_supply = sum == 0)) return;
669
670      // Copying _graph_ref to graph
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), _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      // Copying _graph_ref to graph
709      copyGraph(_graph, _graph_ref)
710        .edgeMap(_cost, cost)
711        .nodeRef(_node_ref)
712        .edgeRef(_edge_ref)
713        .run();
714
715      // Removing non-zero lower bounds
716      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
717        _capacity[_edge_ref[e]] = capacity[e] - lower[e];
718      }
719      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {
720        Supply sum = 0;
721        if (n == s) sum =  flow_value;
722        if (n == t) sum = -flow_value;
723        for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)
724          sum += lower[e];
725        for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)
726          sum -= lower[e];
727        _supply[_node_ref[n]] = sum;
728      }
729      _valid_supply = true;
730    }
731
732    /// \brief Simple constructor (without lower bounds).
733    ///
734    /// Simple constructor (without lower bounds).
735    ///
736    /// \param graph The directed graph the algorithm runs on.
737    /// \param capacity The capacities (upper bounds) of the edges.
738    /// \param cost The cost (length) values of the edges.
739    /// \param s The source node.
740    /// \param t The target node.
741    /// \param flow_value The required amount of flow from node \c s
742    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
743    NetworkSimplex( const Graph &graph,
744                    const CapacityMap &capacity,
745                    const CostMap &cost,
746                    typename Graph::Node s,
747                    typename Graph::Node t,
748                    typename SupplyMap::Value flow_value ) :
749      _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
750      _cost(_graph), _supply(_graph, 0), _flow(_graph),
751      _potential(_graph), _depth(_graph), _parent(_graph),
752      _pred_edge(_graph), _thread(_graph), _forward(_graph),
753      _state(_graph), _red_cost(_graph, _cost, _potential),
754      _flow_result(NULL), _potential_result(NULL),
755      _local_flow(false), _local_potential(false),
756      _node_ref(graph), _edge_ref(graph)
757    {
758      // Copying _graph_ref to graph
759      copyGraph(_graph, _graph_ref)
760        .edgeMap(_capacity, capacity)
761        .edgeMap(_cost, cost)
762        .nodeRef(_node_ref)
763        .edgeRef(_edge_ref)
764        .run();
765      _supply[_node_ref[s]] =  flow_value;
766      _supply[_node_ref[t]] = -flow_value;
767      _valid_supply = true;
768    }
769
770    /// Destructor.
771    ~NetworkSimplex() {
772      if (_local_flow) delete _flow_result;
773      if (_local_potential) delete _potential_result;
774    }
775
776    /// \brief Set the flow map.
777    ///
778    /// Set the flow map.
779    ///
780    /// \return \c (*this)
781    NetworkSimplex& flowMap(FlowMap &map) {
782      if (_local_flow) {
783        delete _flow_result;
784        _local_flow = false;
785      }
786      _flow_result = &map;
787      return *this;
788    }
789
790    /// \brief Set the potential map.
791    ///
792    /// Set the potential map.
793    ///
794    /// \return \c (*this)
795    NetworkSimplex& potentialMap(PotentialMap &map) {
796      if (_local_potential) {
797        delete _potential_result;
798        _local_potential = false;
799      }
800      _potential_result = &map;
801      return *this;
802    }
803
804    /// \name Execution control
805
806    /// @{
807
808    /// \brief Runs the algorithm.
809    ///
810    /// Runs the algorithm.
811    ///
812    /// \param pivot_rule The pivot rule that is used during the
813    /// algorithm.
814    ///
815    /// The available pivot rules:
816    ///
817    /// - FIRST_ELIGIBLE_PIVOT The next eligible edge is selected in
818    /// a wraparound fashion in every iteration
819    /// (\ref FirstEligiblePivotRule).
820    ///
821    /// - BEST_ELIGIBLE_PIVOT The best eligible edge is selected in
822    /// every iteration (\ref BestEligiblePivotRule).
823    ///
824    /// - BLOCK_SEARCH_PIVOT A specified number of edges are examined in
825    /// every iteration in a wraparound fashion and the best eligible
826    /// edge is selected from this block (\ref BlockSearchPivotRule).
827    ///
828    /// - CANDIDATE_LIST_PIVOT In a major iteration a candidate list is
829    /// built from eligible edges in a wraparound fashion and in the
830    /// following minor iterations the best eligible edge is selected
831    /// from this list (\ref CandidateListPivotRule).
832    ///
833    /// - ALTERING_LIST_PIVOT It is a modified version of the
834    /// "Candidate List" pivot rule. It keeps only the several best
835    /// eligible edges from the former candidate list and extends this
836    /// list in every iteration (\ref AlteringListPivotRule).
837    ///
838    /// According to our comprehensive benchmark tests the "Block Search"
839    /// pivot rule proved to be by far the fastest and the most robust
840    /// on various test inputs. Thus it is the default option.
841    ///
842    /// \return \c true if a feasible flow can be found.
843    bool run(PivotRuleEnum pivot_rule = BLOCK_SEARCH_PIVOT) {
844      return init() && start(pivot_rule);
845    }
846
847    /// @}
848
849    /// \name Query Functions
850    /// The results of the algorithm can be obtained using these
851    /// functions.\n
852    /// \ref lemon::NetworkSimplex::run() "run()" must be called before
853    /// using them.
854
855    /// @{
856
857    /// \brief Return a const reference to the edge map storing the
858    /// found flow.
859    ///
860    /// Return a const reference to the edge map storing the found flow.
861    ///
862    /// \pre \ref run() must be called before using this function.
863    const FlowMap& flowMap() const {
864      return *_flow_result;
865    }
866
867    /// \brief Return a const reference to the node map storing the
868    /// found potentials (the dual solution).
869    ///
870    /// Return a const reference to the node map storing the found
871    /// potentials (the dual solution).
872    ///
873    /// \pre \ref run() must be called before using this function.
874    const PotentialMap& potentialMap() const {
875      return *_potential_result;
876    }
877
878    /// \brief Return the flow on the given edge.
879    ///
880    /// Return the flow on the given edge.
881    ///
882    /// \pre \ref run() must be called before using this function.
883    Capacity flow(const typename Graph::Edge& edge) const {
884      return (*_flow_result)[edge];
885    }
886
887    /// \brief Return the potential of the given node.
888    ///
889    /// Return the potential of the given node.
890    ///
891    /// \pre \ref run() must be called before using this function.
892    Cost potential(const typename Graph::Node& node) const {
893      return (*_potential_result)[node];
894    }
895
896    /// \brief Return the total cost of the found flow.
897    ///
898    /// Return the total cost of the found flow. The complexity of the
899    /// function is \f$ O(e) \f$.
900    ///
901    /// \pre \ref run() must be called before using this function.
902    Cost totalCost() const {
903      Cost c = 0;
904      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
905        c += (*_flow_result)[e] * _cost[_edge_ref[e]];
906      return c;
907    }
908
909    /// @}
910
911  private:
912
913    /// \brief Extend the underlying graph and initialize all the
914    /// node and edge maps.
915    bool init() {
916      if (!_valid_supply) return false;
917
918      // Initializing result maps
919      if (!_flow_result) {
920        _flow_result = new FlowMap(_graph_ref);
921        _local_flow = true;
922      }
923      if (!_potential_result) {
924        _potential_result = new PotentialMap(_graph_ref);
925        _local_potential = true;
926      }
927
928      // Initializing the edge vector and the edge maps
929      _edges.reserve(countEdges(_graph));
930      for (EdgeIt e(_graph); e != INVALID; ++e) {
931        _edges.push_back(e);
932        _flow[e] = 0;
933        _state[e] = STATE_LOWER;
934      }
935
936      // Adding an artificial root node to the graph
937      _root = _graph.addNode();
938      _parent[_root] = INVALID;
939      _pred_edge[_root] = INVALID;
940      _depth[_root] = 0;
941      _supply[_root] = 0;
942      _potential[_root] = 0;
943
944      // Adding artificial edges to the graph and initializing the node
945      // maps of the spanning tree data structure
946      Node last = _root;
947      Edge e;
948      Cost max_cost = std::numeric_limits<Cost>::max() / 4;
949      for (NodeIt u(_graph); u != INVALID; ++u) {
950        if (u == _root) continue;
951        _thread[last] = u;
952        last = u;
953        _parent[u] = _root;
954        _depth[u] = 1;
955        if (_supply[u] >= 0) {
956          e = _graph.addEdge(u, _root);
957          _flow[e] = _supply[u];
958          _forward[u] = true;
959          _potential[u] = -max_cost;
960        } else {
961          e = _graph.addEdge(_root, u);
962          _flow[e] = -_supply[u];
963          _forward[u] = false;
964          _potential[u] = max_cost;
965        }
966        _cost[e] = max_cost;
967        _capacity[e] = std::numeric_limits<Capacity>::max();
968        _state[e] = STATE_TREE;
969        _pred_edge[u] = e;
970      }
971      _thread[last] = _root;
972
973      return true;
974    }
975
976    /// Find the join node.
977    inline Node findJoinNode() {
978      Node u = _graph.source(_in_edge);
979      Node v = _graph.target(_in_edge);
980      while (u != v) {
981        if (_depth[u] == _depth[v]) {
982          u = _parent[u];
983          v = _parent[v];
984        }
985        else if (_depth[u] > _depth[v]) u = _parent[u];
986        else v = _parent[v];
987      }
988      return u;
989    }
990
991    /// \brief Find the leaving edge of the cycle.
992    /// \return \c true if the leaving edge is not the same as the
993    /// entering edge.
994    inline bool findLeavingEdge() {
995      // Initializing first and second nodes according to the direction
996      // of the cycle
997      if (_state[_in_edge] == STATE_LOWER) {
998        first  = _graph.source(_in_edge);
999        second = _graph.target(_in_edge);
1000      } else {
1001        first  = _graph.target(_in_edge);
1002        second = _graph.source(_in_edge);
1003      }
1004      delta = _capacity[_in_edge];
1005      bool result = false;
1006      Capacity d;
1007      Edge e;
1008
1009      // Searching the cycle along the path form the first node to the
1010      // root node
1011      for (Node u = first; u != join; u = _parent[u]) {
1012        e = _pred_edge[u];
1013        d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e];
1014        if (d < delta) {
1015          delta = d;
1016          u_out = u;
1017          u_in = first;
1018          v_in = second;
1019          result = true;
1020        }
1021      }
1022      // Searching the cycle along the path form the second node to the
1023      // root node
1024      for (Node u = second; u != join; u = _parent[u]) {
1025        e = _pred_edge[u];
1026        d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e];
1027        if (d <= delta) {
1028          delta = d;
1029          u_out = u;
1030          u_in = second;
1031          v_in = first;
1032          result = true;
1033        }
1034      }
1035      return result;
1036    }
1037
1038    /// Change \c flow and \c state edge maps.
1039    inline void changeFlows(bool change) {
1040      // Augmenting along the cycle
1041      if (delta > 0) {
1042        Capacity val = _state[_in_edge] * delta;
1043        _flow[_in_edge] += val;
1044        for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) {
1045          _flow[_pred_edge[u]] += _forward[u] ? -val : val;
1046        }
1047        for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) {
1048          _flow[_pred_edge[u]] += _forward[u] ? val : -val;
1049        }
1050      }
1051      // Updating the state of the entering and leaving edges
1052      if (change) {
1053        _state[_in_edge] = STATE_TREE;
1054        _state[_pred_edge[u_out]] =
1055          (_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
1056      } else {
1057        _state[_in_edge] = -_state[_in_edge];
1058      }
1059    }
1060
1061    /// Update \c thread and \c parent node maps.
1062    inline void updateThreadParent() {
1063      Node u;
1064      v_out = _parent[u_out];
1065
1066      // Handling the case when join and v_out coincide
1067      bool par_first = false;
1068      if (join == v_out) {
1069        for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
1070        if (u == v_in) {
1071          par_first = true;
1072          while (_thread[u] != u_out) u = _thread[u];
1073          first = u;
1074        }
1075      }
1076
1077      // Finding the last successor of u_in (u) and the node after it
1078      // (right) according to the thread index
1079      for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
1080      right = _thread[u];
1081      if (_thread[v_in] == u_out) {
1082        for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;
1083        if (last == u_out) last = _thread[last];
1084      }
1085      else last = _thread[v_in];
1086
1087      // Updating stem nodes
1088      _thread[v_in] = stem = u_in;
1089      par_stem = v_in;
1090      while (stem != u_out) {
1091        _thread[u] = new_stem = _parent[stem];
1092
1093        // Finding the node just before the stem node (u) according to
1094        // the original thread index
1095        for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
1096        _thread[u] = right;
1097
1098        // Changing the parent node of stem and shifting stem and
1099        // par_stem nodes
1100        _parent[stem] = par_stem;
1101        par_stem = stem;
1102        stem = new_stem;
1103
1104        // Finding the last successor of stem (u) and the node after it
1105        // (right) according to the thread index
1106        for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
1107        right = _thread[u];
1108      }
1109      _parent[u_out] = par_stem;
1110      _thread[u] = last;
1111
1112      if (join == v_out && par_first) {
1113        if (first != v_in) _thread[first] = right;
1114      } else {
1115        for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
1116        _thread[u] = right;
1117      }
1118    }
1119
1120    /// Update \c pred_edge and \c forward node maps.
1121    inline void updatePredEdge() {
1122      Node u = u_out, v;
1123      while (u != u_in) {
1124        v = _parent[u];
1125        _pred_edge[u] = _pred_edge[v];
1126        _forward[u] = !_forward[v];
1127        u = v;
1128      }
1129      _pred_edge[u_in] = _in_edge;
1130      _forward[u_in] = (u_in == _graph.source(_in_edge));
1131    }
1132
1133    /// Update \c depth and \c potential node maps.
1134    inline void updateDepthPotential() {
1135      _depth[u_in] = _depth[v_in] + 1;
1136      _potential[u_in] = _forward[u_in] ?
1137        _potential[v_in] - _cost[_pred_edge[u_in]] :
1138        _potential[v_in] + _cost[_pred_edge[u_in]];
1139
1140      Node u = _thread[u_in], v;
1141      while (true) {
1142        v = _parent[u];
1143        if (v == INVALID) break;
1144        _depth[u] = _depth[v] + 1;
1145        _potential[u] = _forward[u] ?
1146          _potential[v] - _cost[_pred_edge[u]] :
1147          _potential[v] + _cost[_pred_edge[u]];
1148        if (_depth[u] <= _depth[v_in]) break;
1149        u = _thread[u];
1150      }
1151    }
1152
1153    /// Execute the algorithm.
1154    bool start(PivotRuleEnum pivot_rule) {
1155      // Selecting the pivot rule implementation
1156      switch (pivot_rule) {
1157        case FIRST_ELIGIBLE_PIVOT:
1158          return start<FirstEligiblePivotRule>();
1159        case BEST_ELIGIBLE_PIVOT:
1160          return start<BestEligiblePivotRule>();
1161        case BLOCK_SEARCH_PIVOT:
1162          return start<BlockSearchPivotRule>();
1163        case CANDIDATE_LIST_PIVOT:
1164          return start<CandidateListPivotRule>();
1165        case ALTERING_LIST_PIVOT:
1166          return start<AlteringListPivotRule>();
1167      }
1168      return false;
1169    }
1170
1171    template<class PivotRuleImplementation>
1172    bool start() {
1173      PivotRuleImplementation pivot(*this, _edges);
1174
1175      // Executing the network simplex algorithm
1176      while (pivot.findEnteringEdge()) {
1177        join = findJoinNode();
1178        bool change = findLeavingEdge();
1179        changeFlows(change);
1180        if (change) {
1181          updateThreadParent();
1182          updatePredEdge();
1183          updateDepthPotential();
1184        }
1185      }
1186
1187      // Checking if the flow amount equals zero on all the artificial
1188      // edges
1189      for (InEdgeIt e(_graph, _root); e != INVALID; ++e)
1190        if (_flow[e] > 0) return false;
1191      for (OutEdgeIt e(_graph, _root); e != INVALID; ++e)
1192        if (_flow[e] > 0) return false;
1193
1194      // Copying flow values to _flow_result
1195      if (_lower) {
1196        for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
1197          (*_flow_result)[e] = (*_lower)[e] + _flow[_edge_ref[e]];
1198      } else {
1199        for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
1200          (*_flow_result)[e] = _flow[_edge_ref[e]];
1201      }
1202      // Copying potential values to _potential_result
1203      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
1204        (*_potential_result)[n] = _potential[_node_ref[n]];
1205
1206      return true;
1207    }
1208
1209  }; //class NetworkSimplex
1210
1211  ///@}
1212
1213} //namespace lemon
1214
1215#endif //LEMON_NETWORK_SIMPLEX_H
Note: See TracBrowser for help on using the repository browser.