COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/network_simplex.h @ 2533:aea952a1af99

Last change on this file since 2533:aea952a1af99 was 2533:aea952a1af99, checked in by Peter Kovacs, 12 years ago

Bug fixes.

File size: 27.4 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2007
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 The network simplex algorithm for finding a minimum cost
26/// flow.
27
28#include <limits>
29#include <lemon/graph_adaptor.h>
30#include <lemon/graph_utils.h>
31#include <lemon/smart_graph.h>
32
33/// \brief The pivot rule used in the algorithm.
34//#define FIRST_ELIGIBLE_PIVOT
35//#define BEST_ELIGIBLE_PIVOT
36#define EDGE_BLOCK_PIVOT
37//#define CANDIDATE_LIST_PIVOT
38//#define SORTED_LIST_PIVOT
39
40//#define _DEBUG_ITER_
41
42
43/// \brief State constant for edges at their lower bounds.
44#define LOWER   1
45/// \brief State constant for edges in the spanning tree.
46#define TREE    0
47/// \brief State constant for edges at their upper bounds.
48#define UPPER   -1
49
50#ifdef EDGE_BLOCK_PIVOT
51  #include <cmath>
52  /// \brief Lower bound for the size of blocks.
53  #define MIN_BLOCK_SIZE        10
54#endif
55
56#ifdef CANDIDATE_LIST_PIVOT
57  #include <vector>
58  #define LIST_LENGTH_DIV           50
59  #define MINOR_LIMIT_DIV           200
60#endif
61
62#ifdef SORTED_LIST_PIVOT
63  #include <vector>
64  #include <algorithm>
65  #define LIST_LENGTH_DIV       100
66  #define LOWER_DIV             4
67#endif
68
69namespace lemon {
70
71  /// \addtogroup min_cost_flow
72  /// @{
73
74  /// \brief Implementation of the network simplex algorithm for
75  /// finding a minimum cost flow.
76  ///
77  /// \ref lemon::NetworkSimplex "NetworkSimplex" implements the
78  /// network simplex algorithm for finding a minimum cost flow.
79  ///
80  /// \param Graph The directed graph type the algorithm runs on.
81  /// \param LowerMap The type of the lower bound map.
82  /// \param CapacityMap The type of the capacity (upper bound) map.
83  /// \param CostMap The type of the cost (length) map.
84  /// \param SupplyMap The type of the supply map.
85  ///
86  /// \warning
87  /// - Edge capacities and costs should be nonnegative integers.
88  ///   However \c CostMap::Value should be signed type.
89  /// - Supply values should be signed integers.
90  /// - \c LowerMap::Value must be convertible to
91  ///   \c CapacityMap::Value and \c CapacityMap::Value must be
92  ///   convertible to \c SupplyMap::Value.
93  ///
94  /// \author Peter Kovacs
95
96  template < typename Graph,
97             typename LowerMap = typename Graph::template EdgeMap<int>,
98             typename CapacityMap = LowerMap,
99             typename CostMap = typename Graph::template EdgeMap<int>,
100             typename SupplyMap = typename Graph::template NodeMap
101                                  <typename CapacityMap::Value> >
102  class NetworkSimplex
103  {
104    typedef typename LowerMap::Value Lower;
105    typedef typename CapacityMap::Value Capacity;
106    typedef typename CostMap::Value Cost;
107    typedef typename SupplyMap::Value Supply;
108
109    typedef SmartGraph SGraph;
110    typedef typename SGraph::Node Node;
111    typedef typename SGraph::NodeIt NodeIt;
112    typedef typename SGraph::Edge Edge;
113    typedef typename SGraph::EdgeIt EdgeIt;
114    typedef typename SGraph::InEdgeIt InEdgeIt;
115    typedef typename SGraph::OutEdgeIt OutEdgeIt;
116
117    typedef typename SGraph::template EdgeMap<Lower> SLowerMap;
118    typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap;
119    typedef typename SGraph::template EdgeMap<Cost> SCostMap;
120    typedef typename SGraph::template NodeMap<Supply> SSupplyMap;
121    typedef typename SGraph::template NodeMap<Cost> SPotentialMap;
122
123    typedef typename SGraph::template NodeMap<int> IntNodeMap;
124    typedef typename SGraph::template NodeMap<bool> BoolNodeMap;
125    typedef typename SGraph::template NodeMap<Node> NodeNodeMap;
126    typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap;
127    typedef typename SGraph::template EdgeMap<int> IntEdgeMap;
128
129    typedef typename Graph::template NodeMap<Node> NodeRefMap;
130    typedef typename Graph::template EdgeMap<Edge> EdgeRefMap;
131
132  public:
133
134    /// \brief The type of the flow map.
135    typedef typename Graph::template EdgeMap<Capacity> FlowMap;
136    /// \brief The type of the potential map.
137    typedef typename Graph::template NodeMap<Cost> PotentialMap;
138
139  protected:
140
141    /// \brief Map adaptor class for handling reduced edge costs.
142    class ReducedCostMap : public MapBase<Edge, Cost>
143    {
144    private:
145
146      const SGraph &gr;
147      const SCostMap &cost_map;
148      const SPotentialMap &pot_map;
149
150    public:
151
152      ReducedCostMap( const SGraph &_gr,
153                      const SCostMap &_cm,
154                      const SPotentialMap &_pm ) :
155        gr(_gr), cost_map(_cm), pot_map(_pm) {}
156
157      Cost operator[](const Edge &e) const {
158        return cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];
159      }
160
161    }; //class ReducedCostMap
162
163  protected:
164
165    /// \brief The directed graph the algorithm runs on.
166    SGraph graph;
167    /// \brief The original graph.
168    const Graph &graph_ref;
169    /// \brief The original lower bound map.
170    const LowerMap *lower;
171    /// \brief The capacity map.
172    SCapacityMap capacity;
173    /// \brief The cost map.
174    SCostMap cost;
175    /// \brief The supply map.
176    SSupplyMap supply;
177    /// \brief The reduced cost map.
178    ReducedCostMap red_cost;
179    /// \brief The sum of supply values equals zero.
180    bool valid_supply;
181
182    /// \brief The edge map of the current flow.
183    SCapacityMap flow;
184    /// \brief The edge map of the found flow on the original graph.
185    FlowMap flow_result;
186    /// \brief The potential node map.
187    SPotentialMap potential;
188    /// \brief The potential node map on the original graph.
189    PotentialMap potential_result;
190
191    /// \brief Node reference for the original graph.
192    NodeRefMap node_ref;
193    /// \brief Edge reference for the original graph.
194    EdgeRefMap edge_ref;
195
196    /// \brief The depth node map of the spanning tree structure.
197    IntNodeMap depth;
198    /// \brief The parent node map of the spanning tree structure.
199    NodeNodeMap parent;
200    /// \brief The pred_edge node map of the spanning tree structure.
201    EdgeNodeMap pred_edge;
202    /// \brief The thread node map of the spanning tree structure.
203    NodeNodeMap thread;
204    /// \brief The forward node map of the spanning tree structure.
205    BoolNodeMap forward;
206    /// \brief The state edge map.
207    IntEdgeMap state;
208
209
210#ifdef EDGE_BLOCK_PIVOT
211    /// \brief The size of blocks for the "Edge Block" pivot rule.
212    int block_size;
213#endif
214#ifdef CANDIDATE_LIST_PIVOT
215    /// \brief The list of candidate edges for the "Candidate List"
216    /// pivot rule.
217    std::vector<Edge> candidates;
218    /// \brief The maximum length of the edge list for the
219    /// "Candidate List" pivot rule.
220    int list_length;
221    /// \brief The maximum number of minor iterations between two major
222    /// itarations.
223    int minor_limit;
224    /// \brief The number of minor iterations.
225    int minor_count;
226#endif
227#ifdef SORTED_LIST_PIVOT
228    /// \brief The list of candidate edges for the "Sorted List"
229    /// pivot rule.
230    std::vector<Edge> candidates;
231    /// \brief The maximum length of the edge list for the
232    /// "Sorted List" pivot rule.
233    int list_length;
234    int list_index;
235#endif
236
237    // Root node of the starting spanning tree.
238    Node root;
239    // The entering edge of the current pivot iteration.
240    Edge in_edge;
241    // Temporary nodes used in the current pivot iteration.
242    Node join, u_in, v_in, u_out, v_out;
243    Node right, first, second, last;
244    Node stem, par_stem, new_stem;
245    // The maximum augment amount along the cycle in the current pivot
246    // iteration.
247    Capacity delta;
248
249  public :
250
251    /// \brief General constructor of the class (with lower bounds).
252    ///
253    /// General constructor of the class (with lower bounds).
254    ///
255    /// \param _graph The directed graph the algorithm runs on.
256    /// \param _lower The lower bounds of the edges.
257    /// \param _capacity The capacities (upper bounds) of the edges.
258    /// \param _cost The cost (length) values of the edges.
259    /// \param _supply The supply values of the nodes (signed).
260    NetworkSimplex( const Graph &_graph,
261                    const LowerMap &_lower,
262                    const CapacityMap &_capacity,
263                    const CostMap &_cost,
264                    const SupplyMap &_supply ) :
265      graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),
266      supply(graph), flow(graph), flow_result(_graph), potential(graph),
267      potential_result(_graph), depth(graph), parent(graph),
268      pred_edge(graph), thread(graph), forward(graph), state(graph),
269      node_ref(graph_ref), edge_ref(graph_ref),
270      red_cost(graph, cost, potential)
271    {
272      // Checking the sum of supply values
273      Supply sum = 0;
274      for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
275        sum += _supply[n];
276      if (!(valid_supply = sum == 0)) return;
277
278      // Copying graph_ref to graph
279      graph.reserveNode(countNodes(graph_ref) + 1);
280      graph.reserveEdge(countEdges(graph_ref) + countNodes(graph_ref));
281      copyGraph(graph, graph_ref)
282        .edgeMap(cost, _cost)
283        .nodeRef(node_ref)
284        .edgeRef(edge_ref)
285        .run();
286
287      // Removing nonzero lower bounds
288      for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {
289        capacity[edge_ref[e]] = _capacity[e] - _lower[e];
290      }
291      for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {
292        Supply s = _supply[n];
293        for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
294          s += _lower[e];
295        for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
296          s -= _lower[e];
297        supply[node_ref[n]] = s;
298      }
299    }
300
301    /// \brief General constructor of the class (without lower bounds).
302    ///
303    /// General constructor of the class (without lower bounds).
304    ///
305    /// \param _graph The directed graph the algorithm runs on.
306    /// \param _capacity The capacities (upper bounds) of the edges.
307    /// \param _cost The cost (length) values of the edges.
308    /// \param _supply The supply values of the nodes (signed).
309    NetworkSimplex( const Graph &_graph,
310                    const CapacityMap &_capacity,
311                    const CostMap &_cost,
312                    const SupplyMap &_supply ) :
313      graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),
314      supply(graph), flow(graph), flow_result(_graph), potential(graph),
315      potential_result(_graph), depth(graph), parent(graph),
316      pred_edge(graph), thread(graph), forward(graph), state(graph),
317      node_ref(graph_ref), edge_ref(graph_ref),
318      red_cost(graph, cost, potential)
319    {
320      // Checking the sum of supply values
321      Supply sum = 0;
322      for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
323        sum += _supply[n];
324      if (!(valid_supply = sum == 0)) return;
325
326      // Copying graph_ref to graph
327      copyGraph(graph, graph_ref)
328        .edgeMap(capacity, _capacity)
329        .edgeMap(cost, _cost)
330        .nodeMap(supply, _supply)
331        .nodeRef(node_ref)
332        .edgeRef(edge_ref)
333        .run();
334    }
335
336    /// \brief Simple constructor of the class (with lower bounds).
337    ///
338    /// Simple constructor of the class (with lower bounds).
339    ///
340    /// \param _graph The directed graph the algorithm runs on.
341    /// \param _lower The lower bounds of the edges.
342    /// \param _capacity The capacities (upper bounds) of the edges.
343    /// \param _cost The cost (length) values of the edges.
344    /// \param _s The source node.
345    /// \param _t The target node.
346    /// \param _flow_value The required amount of flow from node \c _s
347    /// to node \c _t (i.e. the supply of \c _s and the demand of
348    /// \c _t).
349    NetworkSimplex( const Graph &_graph,
350                    const LowerMap &_lower,
351                    const CapacityMap &_capacity,
352                    const CostMap &_cost,
353                    typename Graph::Node _s,
354                    typename Graph::Node _t,
355                    typename SupplyMap::Value _flow_value ) :
356      graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),
357      supply(graph), flow(graph), flow_result(_graph), potential(graph),
358      potential_result(_graph), depth(graph), parent(graph),
359      pred_edge(graph), thread(graph), forward(graph), state(graph),
360      node_ref(graph_ref), edge_ref(graph_ref),
361      red_cost(graph, cost, potential)
362    {
363      // Copying graph_ref to graph
364      copyGraph(graph, graph_ref)
365        .edgeMap(cost, _cost)
366        .nodeRef(node_ref)
367        .edgeRef(edge_ref)
368        .run();
369
370      // Removing nonzero lower bounds
371      for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {
372        capacity[edge_ref[e]] = _capacity[e] - _lower[e];
373      }
374      for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {
375        Supply s = 0;
376        if (n == _s) s =  _flow_value;
377        if (n == _t) s = -_flow_value;
378        for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
379          s += _lower[e];
380        for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
381          s -= _lower[e];
382        supply[node_ref[n]] = s;
383      }
384      valid_supply = true;
385    }
386
387    /// \brief Simple constructor of the class (without lower bounds).
388    ///
389    /// Simple constructor of the class (without lower bounds).
390    ///
391    /// \param _graph The directed graph the algorithm runs on.
392    /// \param _capacity The capacities (upper bounds) of the edges.
393    /// \param _cost The cost (length) values of the edges.
394    /// \param _s The source node.
395    /// \param _t The target node.
396    /// \param _flow_value The required amount of flow from node \c _s
397    /// to node \c _t (i.e. the supply of \c _s and the demand of
398    /// \c _t).
399    NetworkSimplex( const Graph &_graph,
400                    const CapacityMap &_capacity,
401                    const CostMap &_cost,
402                    typename Graph::Node _s,
403                    typename Graph::Node _t,
404                    typename SupplyMap::Value _flow_value ) :
405      graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),
406      supply(graph, 0), flow(graph), flow_result(_graph), potential(graph),
407      potential_result(_graph), depth(graph), parent(graph),
408      pred_edge(graph), thread(graph), forward(graph), state(graph),
409      node_ref(graph_ref), edge_ref(graph_ref),
410      red_cost(graph, cost, potential)
411    {
412      // Copying graph_ref to graph
413      copyGraph(graph, graph_ref)
414        .edgeMap(capacity, _capacity)
415        .edgeMap(cost, _cost)
416        .nodeRef(node_ref)
417        .edgeRef(edge_ref)
418        .run();
419      supply[node_ref[_s]] =  _flow_value;
420      supply[node_ref[_t]] = -_flow_value;
421      valid_supply = true;
422    }
423
424    /// \brief Returns a const reference to the flow map.
425    ///
426    /// Returns a const reference to the flow map.
427    ///
428    /// \pre \ref run() must be called before using this function.
429    const FlowMap& flowMap() const {
430      return flow_result;
431    }
432
433    /// \brief Returns a const reference to the potential map (the dual
434    /// solution).
435    ///
436    /// Returns a const reference to the potential map (the dual
437    /// solution).
438    ///
439    /// \pre \ref run() must be called before using this function.
440    const PotentialMap& potentialMap() const {
441      return potential_result;
442    }
443
444    /// \brief Returns the total cost of the found flow.
445    ///
446    /// Returns the total cost of the found flow. The complexity of the
447    /// function is \f$ O(e) \f$.
448    ///
449    /// \pre \ref run() must be called before using this function.
450    Cost totalCost() const {
451      Cost c = 0;
452      for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
453        c += flow_result[e] * cost[edge_ref[e]];
454      return c;
455    }
456
457    /// \brief Runs the algorithm.
458    ///
459    /// Runs the algorithm.
460    ///
461    /// \return \c true if a feasible flow can be found.
462    bool run() {
463      return init() && start();
464    }
465
466  protected:
467
468    /// \brief Extends the underlaying graph and initializes all the
469    /// node and edge maps.
470    bool init() {
471      if (!valid_supply) return false;
472
473      // Initializing state and flow maps
474      for (EdgeIt e(graph); e != INVALID; ++e) {
475        flow[e] = 0;
476        state[e] = LOWER;
477      }
478
479      // Adding an artificial root node to the graph
480      root = graph.addNode();
481      parent[root] = INVALID;
482      pred_edge[root] = INVALID;
483      depth[root] = 0;
484      supply[root] = 0;
485      potential[root] = 0;
486
487      // Adding artificial edges to the graph and initializing the node
488      // maps of the spanning tree data structure
489      Supply sum = 0;
490      Node last = root;
491      Edge e;
492      Cost max_cost = std::numeric_limits<Cost>::max() / 4;
493      for (NodeIt u(graph); u != INVALID; ++u) {
494        if (u == root) continue;
495        thread[last] = u;
496        last = u;
497        parent[u] = root;
498        depth[u] = 1;
499        sum += supply[u];
500        if (supply[u] >= 0) {
501          e = graph.addEdge(u, root);
502          flow[e] = supply[u];
503          forward[u] = true;
504          potential[u] = max_cost;
505        } else {
506          e = graph.addEdge(root, u);
507          flow[e] = -supply[u];
508          forward[u] = false;
509          potential[u] = -max_cost;
510        }
511        cost[e] = max_cost;
512        capacity[e] = std::numeric_limits<Capacity>::max();
513        state[e] = TREE;
514        pred_edge[u] = e;
515      }
516      thread[last] = root;
517
518#ifdef EDGE_BLOCK_PIVOT
519      // Initializing block_size for the edge block pivot rule
520      int edge_num = countEdges(graph);
521      block_size = 2 * int(sqrt(countEdges(graph)));
522      if (block_size < MIN_BLOCK_SIZE) block_size = MIN_BLOCK_SIZE;
523//      block_size = edge_num >= BLOCK_NUM * MIN_BLOCK_SIZE ?
524//                   edge_num / BLOCK_NUM : MIN_BLOCK_SIZE;
525#endif
526#ifdef CANDIDATE_LIST_PIVOT
527      int edge_num = countEdges(graph);
528      minor_count = 0;
529      list_length = edge_num / LIST_LENGTH_DIV;
530      minor_limit = edge_num / MINOR_LIMIT_DIV;
531#endif
532#ifdef SORTED_LIST_PIVOT
533      int edge_num = countEdges(graph);
534      list_index = 0;
535      list_length = edge_num / LIST_LENGTH_DIV;
536#endif
537
538      return sum == 0;
539    }
540
541#ifdef FIRST_ELIGIBLE_PIVOT
542    /// \brief Finds entering edge according to the "First Eligible"
543    /// pivot rule.
544    bool findEnteringEdge(EdgeIt &next_edge) {
545      for (EdgeIt e = next_edge; e != INVALID; ++e) {
546        if (state[e] * red_cost[e] < 0) {
547          in_edge = e;
548          next_edge = ++e;
549          return true;
550        }
551      }
552      for (EdgeIt e(graph); e != next_edge; ++e) {
553        if (state[e] * red_cost[e] < 0) {
554          in_edge = e;
555          next_edge = ++e;
556          return true;
557        }
558      }
559      return false;
560    }
561#endif
562
563#ifdef BEST_ELIGIBLE_PIVOT
564    /// \brief Finds entering edge according to the "Best Eligible"
565    /// pivot rule.
566    bool findEnteringEdge() {
567      Cost min = 0;
568      for (EdgeIt e(graph); e != INVALID; ++e) {
569        if (state[e] * red_cost[e] < min) {
570          min = state[e] * red_cost[e];
571          in_edge = e;
572        }
573      }
574      return min < 0;
575    }
576#endif
577
578#ifdef EDGE_BLOCK_PIVOT
579    /// \brief Finds entering edge according to the "Edge Block"
580    /// pivot rule.
581    bool findEnteringEdge(EdgeIt &next_edge) {
582      // Performing edge block selection
583      Cost curr, min = 0;
584      EdgeIt min_edge(graph);
585      int cnt = 0;
586      for (EdgeIt e = next_edge; e != INVALID; ++e) {
587        if ((curr = state[e] * red_cost[e]) < min) {
588          min = curr;
589          min_edge = e;
590        }
591        if (++cnt == block_size) {
592          if (min < 0) break;
593          cnt = 0;
594        }
595      }
596      if (!(min < 0)) {
597        for (EdgeIt e(graph); e != next_edge; ++e) {
598          if ((curr = state[e] * red_cost[e]) < min) {
599            min = curr;
600            min_edge = e;
601          }
602          if (++cnt == block_size) {
603            if (min < 0) break;
604            cnt = 0;
605          }
606        }
607      }
608      in_edge = min_edge;
609      if ((next_edge = ++min_edge) == INVALID)
610        next_edge = EdgeIt(graph);
611      return min < 0;
612    }
613#endif
614
615#ifdef CANDIDATE_LIST_PIVOT
616    /// \brief Finds entering edge according to the "Candidate List"
617    /// pivot rule.
618    bool findEnteringEdge() {
619      typedef typename std::vector<Edge>::iterator ListIt;
620
621      if (minor_count >= minor_limit || candidates.size() == 0) {
622        // Major iteration
623        candidates.clear();
624        for (EdgeIt e(graph); e != INVALID; ++e) {
625          if (state[e] * red_cost[e] < 0) {
626            candidates.push_back(e);
627            if (candidates.size() == list_length) break;
628          }
629        }
630        if (candidates.size() == 0) return false;
631      }
632
633      // Minor iteration
634      ++minor_count;
635      Cost min = 0;
636      Edge e;
637      for (int i = 0; i < candidates.size(); ++i) {
638        e = candidates[i];
639        if (state[e] * red_cost[e] < min) {
640          min = state[e] * red_cost[e];
641          in_edge = e;
642        }
643      }
644      return true;
645    }
646#endif
647
648#ifdef SORTED_LIST_PIVOT
649    /// \brief Functor class to compare edges during sort of the
650    /// candidate list.
651    class SortFunc
652    {
653    private:
654      const IntEdgeMap &st;
655      const ReducedCostMap &rc;
656    public:
657      SortFunc(const IntEdgeMap &_st, const ReducedCostMap &_rc) :
658        st(_st), rc(_rc) {}
659      bool operator()(const Edge &e1, const Edge &e2) {
660        return st[e1] * rc[e1] < st[e2] * rc[e2];
661      }
662    };
663
664    /// \brief Finds entering edge according to the "Sorted List"
665    /// pivot rule.
666    bool findEnteringEdge() {
667      static SortFunc sort_func(state, red_cost);
668
669      // Minor iteration
670      while (list_index < candidates.size()) {
671        in_edge = candidates[list_index++];
672        if (state[in_edge] * red_cost[in_edge] < 0) return true;
673      }
674
675      // Major iteration
676      candidates.clear();
677      Cost curr, min = 0;
678      for (EdgeIt e(graph); e != INVALID; ++e) {
679        if ((curr = state[e] * red_cost[e]) < min / LOWER_DIV) {
680          candidates.push_back(e);
681          if (curr < min) min = curr;
682          if (candidates.size() == list_length) break;
683        }
684      }
685      if (candidates.size() == 0) return false;
686      sort(candidates.begin(), candidates.end(), sort_func);
687      in_edge = candidates[0];
688      list_index = 1;
689      return true;
690    }
691#endif
692
693    /// \brief Finds the join node.
694    Node findJoinNode() {
695      // Finding the join node
696      Node u = graph.source(in_edge);
697      Node v = graph.target(in_edge);
698      while (u != v) {
699        if (depth[u] == depth[v]) {
700          u = parent[u];
701          v = parent[v];
702        }
703        else if (depth[u] > depth[v]) u = parent[u];
704        else v = parent[v];
705      }
706      return u;
707    }
708
709    /// \brief Finds the leaving edge of the cycle. Returns \c true if
710    /// the leaving edge is not the same as the entering edge.
711    bool findLeavingEdge() {
712      // Initializing first and second nodes according to the direction
713      // of the cycle
714      if (state[in_edge] == LOWER) {
715        first = graph.source(in_edge);
716        second  = graph.target(in_edge);
717      } else {
718        first = graph.target(in_edge);
719        second  = graph.source(in_edge);
720      }
721      delta = capacity[in_edge];
722      bool result = false;
723      Capacity d;
724      Edge e;
725
726      // Searching the cycle along the path form the first node to the
727      // root node
728      for (Node u = first; u != join; u = parent[u]) {
729        e = pred_edge[u];
730        d = forward[u] ? flow[e] : capacity[e] - flow[e];
731        if (d < delta) {
732          delta = d;
733          u_out = u;
734          u_in = first;
735          v_in = second;
736          result = true;
737        }
738      }
739      // Searching the cycle along the path form the second node to the
740      // root node
741      for (Node u = second; u != join; u = parent[u]) {
742        e = pred_edge[u];
743        d = forward[u] ? capacity[e] - flow[e] : flow[e];
744        if (d <= delta) {
745          delta = d;
746          u_out = u;
747          u_in = second;
748          v_in = first;
749          result = true;
750        }
751      }
752      return result;
753    }
754
755    /// \brief Changes flow and state edge maps.
756    void changeFlows(bool change) {
757      // Augmenting along the cycle
758      if (delta > 0) {
759        Capacity val = state[in_edge] * delta;
760        flow[in_edge] += val;
761        for (Node u = graph.source(in_edge); u != join; u = parent[u]) {
762          flow[pred_edge[u]] += forward[u] ? -val : val;
763        }
764        for (Node u = graph.target(in_edge); u != join; u = parent[u]) {
765          flow[pred_edge[u]] += forward[u] ? val : -val;
766        }
767      }
768      // Updating the state of the entering and leaving edges
769      if (change) {
770        state[in_edge] = TREE;
771        state[pred_edge[u_out]] =
772          (flow[pred_edge[u_out]] == 0) ? LOWER : UPPER;
773      } else {
774        state[in_edge] = -state[in_edge];
775      }
776    }
777
778    /// \brief Updates thread and parent node maps.
779    void updateThreadParent() {
780      Node u;
781      v_out = parent[u_out];
782
783      // Handling the case when join and v_out coincide
784      bool par_first = false;
785      if (join == v_out) {
786        for (u = join; u != u_in && u != v_in; u = thread[u]) ;
787        if (u == v_in) {
788          par_first = true;
789          while (thread[u] != u_out) u = thread[u];
790          first = u;
791        }
792      }
793
794      // Finding the last successor of u_in (u) and the node after it
795      // (right) according to the thread index
796      for (u = u_in; depth[thread[u]] > depth[u_in]; u = thread[u]) ;
797      right = thread[u];
798      if (thread[v_in] == u_out) {
799        for (last = u; depth[last] > depth[u_out]; last = thread[last]) ;
800        if (last == u_out) last = thread[last];
801      }
802      else last = thread[v_in];
803
804      // Updating stem nodes
805      thread[v_in] = stem = u_in;
806      par_stem = v_in;
807      while (stem != u_out) {
808        thread[u] = new_stem = parent[stem];
809
810        // Finding the node just before the stem node (u) according to
811        // the original thread index
812        for (u = new_stem; thread[u] != stem; u = thread[u]) ;
813        thread[u] = right;
814
815        // Changing the parent node of stem and shifting stem and
816        // par_stem nodes
817        parent[stem] = par_stem;
818        par_stem = stem;
819        stem = new_stem;
820
821        // Finding the last successor of stem (u) and the node after it
822        // (right) according to the thread index
823        for (u = stem; depth[thread[u]] > depth[stem]; u = thread[u]) ;
824        right = thread[u];
825      }
826      parent[u_out] = par_stem;
827      thread[u] = last;
828
829      if (join == v_out && par_first) {
830        if (first != v_in) thread[first] = right;
831      } else {
832        for (u = v_out; thread[u] != u_out; u = thread[u]) ;
833        thread[u] = right;
834      }
835    }
836
837    /// \brief Updates pred_edge and forward node maps.
838    void updatePredEdge() {
839      Node u = u_out, v;
840      while (u != u_in) {
841        v = parent[u];
842        pred_edge[u] = pred_edge[v];
843        forward[u] = !forward[v];
844        u = v;
845      }
846      pred_edge[u_in] = in_edge;
847      forward[u_in] = (u_in == graph.source(in_edge));
848    }
849
850    /// \brief Updates depth and potential node maps.
851    void updateDepthPotential() {
852      depth[u_in] = depth[v_in] + 1;
853      potential[u_in] = forward[u_in] ?
854        potential[v_in] + cost[pred_edge[u_in]] :
855        potential[v_in] - cost[pred_edge[u_in]];
856
857      Node u = thread[u_in], v;
858      while (true) {
859        v = parent[u];
860        if (v == INVALID) break;
861        depth[u] = depth[v] + 1;
862        potential[u] = forward[u] ?
863          potential[v] + cost[pred_edge[u]] :
864          potential[v] - cost[pred_edge[u]];
865        if (depth[u] <= depth[v_in]) break;
866        u = thread[u];
867      }
868    }
869
870    /// \brief Executes the algorithm.
871    bool start() {
872      // Processing pivots
873#ifdef _DEBUG_ITER_
874      int iter_num = 0;
875#endif
876#if defined(FIRST_ELIGIBLE_PIVOT) || defined(EDGE_BLOCK_PIVOT)
877      EdgeIt next_edge(graph);
878      while (findEnteringEdge(next_edge))
879#else
880      while (findEnteringEdge())
881#endif
882      {
883        join = findJoinNode();
884        bool change = findLeavingEdge();
885        changeFlows(change);
886        if (change) {
887          updateThreadParent();
888          updatePredEdge();
889          updateDepthPotential();
890        }
891#ifdef _DEBUG_ITER_
892        ++iter_num;
893#endif
894      }
895
896#ifdef _DEBUG_ITER_
897      std::cout << "Network Simplex algorithm finished. " << iter_num
898                << " pivot iterations performed." << std::endl;
899#endif
900
901      // Checking if the flow amount equals zero on all the
902      // artificial edges
903      for (InEdgeIt e(graph, root); e != INVALID; ++e)
904        if (flow[e] > 0) return false;
905      for (OutEdgeIt e(graph, root); e != INVALID; ++e)
906        if (flow[e] > 0) return false;
907
908      // Copying flow values to flow_result
909      if (lower) {
910        for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
911          flow_result[e] = (*lower)[e] + flow[edge_ref[e]];
912      } else {
913        for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
914          flow_result[e] = flow[edge_ref[e]];
915      }
916      // Copying potential values to potential_result
917      for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
918        potential_result[n] = potential[node_ref[n]];
919
920      return true;
921    }
922
923  }; //class NetworkSimplex
924
925  ///@}
926
927} //namespace lemon
928
929#endif //LEMON_NETWORK_SIMPLEX_H
Note: See TracBrowser for help on using the repository browser.