COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/cycle_canceling.h @ 2573:a9758ea1f01c

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

Improvements in CycleCanceling?.

Main changes:

  • Use function parameter instead of #define commands to select negative

cycle detection method.

  • Change the name of private members to start with "_".
  • Change the name of function parameters not to start with "_".
  • Remove unnecessary documentation for private members.
  • Doc improvements.
File size: 15.3 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2008
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
12 *
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_CYCLE_CANCELING_H
20#define LEMON_CYCLE_CANCELING_H
21
22/// \ingroup min_cost_flow
23///
24/// \file
25/// \brief Cycle-canceling algorithm for finding a minimum cost flow.
26
27#include <vector>
28#include <lemon/graph_adaptor.h>
29#include <lemon/path.h>
30
31#include <lemon/circulation.h>
32#include <lemon/bellman_ford.h>
33#include <lemon/min_mean_cycle.h>
34
35namespace lemon {
36
37  /// \addtogroup min_cost_flow
38  /// @{
39
40  /// \brief Implementation of a cycle-canceling algorithm for
41  /// finding a minimum cost flow.
42  ///
43  /// \ref CycleCanceling implements a cycle-canceling algorithm for
44  /// finding a minimum cost flow.
45  ///
46  /// \tparam Graph The directed graph type the algorithm runs on.
47  /// \tparam LowerMap The type of the lower bound map.
48  /// \tparam CapacityMap The type of the capacity (upper bound) map.
49  /// \tparam CostMap The type of the cost (length) map.
50  /// \tparam SupplyMap The type of the supply map.
51  ///
52  /// \warning
53  /// - Edge capacities and costs should be \e non-negative \e integers.
54  /// - Supply values should be \e signed \e integers.
55  /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.
56  /// - \c CapacityMap::Value and \c SupplyMap::Value must be
57  ///   convertible to each other.
58  /// - All value types must be convertible to \c CostMap::Value, which
59  ///   must be signed type.
60  ///
61  /// \note By default the \ref BellmanFord "Bellman-Ford" algorithm is
62  /// used for negative cycle detection with limited iteration number.
63  /// However \ref CycleCanceling also provides the "Minimum Mean
64  /// Cycle-Canceling" algorithm, which is \e strongly \e polynomial,
65  /// but rather slower in practice.
66  /// To use this version of the algorithm, call \ref run() with \c true
67  /// parameter.
68  ///
69  /// \author Peter Kovacs
70
71  template < typename Graph,
72             typename LowerMap = typename Graph::template EdgeMap<int>,
73             typename CapacityMap = typename Graph::template EdgeMap<int>,
74             typename CostMap = typename Graph::template EdgeMap<int>,
75             typename SupplyMap = typename Graph::template NodeMap<int> >
76  class CycleCanceling
77  {
78    GRAPH_TYPEDEFS(typename Graph);
79
80    typedef typename CapacityMap::Value Capacity;
81    typedef typename CostMap::Value Cost;
82    typedef typename SupplyMap::Value Supply;
83    typedef typename Graph::template EdgeMap<Capacity> CapacityEdgeMap;
84    typedef typename Graph::template NodeMap<Supply> SupplyNodeMap;
85
86    typedef ResGraphAdaptor< const Graph, Capacity,
87                             CapacityEdgeMap, CapacityEdgeMap > ResGraph;
88    typedef typename ResGraph::Node ResNode;
89    typedef typename ResGraph::NodeIt ResNodeIt;
90    typedef typename ResGraph::Edge ResEdge;
91    typedef typename ResGraph::EdgeIt ResEdgeIt;
92
93  public:
94
95    /// The type of the flow map.
96    typedef typename Graph::template EdgeMap<Capacity> FlowMap;
97
98  private:
99
100    /// \brief Map adaptor class for handling residual edge costs.
101    ///
102    /// \ref ResidualCostMap is a map adaptor class for handling
103    /// residual edge costs.
104    class ResidualCostMap : public MapBase<ResEdge, Cost>
105    {
106    private:
107
108      const CostMap &_cost_map;
109
110    public:
111
112      ///\e
113      ResidualCostMap(const CostMap &cost_map) : _cost_map(cost_map) {}
114
115      ///\e
116      Cost operator[](const ResEdge &e) const {
117        return ResGraph::forward(e) ? _cost_map[e] : -_cost_map[e];
118      }
119
120    }; //class ResidualCostMap
121
122  private:
123
124    // The maximum number of iterations for the first execution of the
125    // Bellman-Ford algorithm. It should be at least 2.
126    static const int BF_FIRST_LIMIT = 2;
127    // The iteration limit for the Bellman-Ford algorithm is multiplied
128    // by BF_ALPHA in every round.
129    static const double BF_ALPHA = 1.5;
130
131  private:
132
133    // The directed graph the algorithm runs on
134    const Graph &_graph;
135    // The original lower bound map
136    const LowerMap *_lower;
137    // The modified capacity map
138    CapacityEdgeMap _capacity;
139    // The original cost map
140    const CostMap &_cost;
141    // The modified supply map
142    SupplyNodeMap _supply;
143    bool _valid_supply;
144
145    // Edge map of the current flow
146    FlowMap _flow;
147
148    // The residual graph
149    ResGraph _res_graph;
150    // The residual cost map
151    ResidualCostMap _res_cost;
152
153  public:
154
155    /// \brief General constructor of the class (with lower bounds).
156    ///
157    /// General constructor of the class (with lower bounds).
158    ///
159    /// \param graph The directed graph the algorithm runs on.
160    /// \param lower The lower bounds of the edges.
161    /// \param capacity The capacities (upper bounds) of the edges.
162    /// \param cost The cost (length) values of the edges.
163    /// \param supply The supply values of the nodes (signed).
164    CycleCanceling( const Graph &graph,
165                    const LowerMap &lower,
166                    const CapacityMap &capacity,
167                    const CostMap &cost,
168                    const SupplyMap &supply ) :
169      _graph(graph), _lower(&lower), _capacity(graph), _cost(cost),
170      _supply(graph), _flow(graph, 0),
171      _res_graph(graph, _capacity, _flow), _res_cost(_cost)
172    {
173      // Removing non-zero lower bounds
174      _capacity = subMap(capacity, lower);
175      Supply sum = 0;
176      for (NodeIt n(_graph); n != INVALID; ++n) {
177        Supply s = supply[n];
178        for (InEdgeIt e(_graph, n); e != INVALID; ++e)
179          s += lower[e];
180        for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
181          s -= lower[e];
182        sum += (_supply[n] = s);
183      }
184      _valid_supply = sum == 0;
185    }
186
187    /// \brief General constructor of the class (without lower bounds).
188    ///
189    /// General constructor of the class (without lower bounds).
190    ///
191    /// \param graph The directed graph the algorithm runs on.
192    /// \param capacity The capacities (upper bounds) of the edges.
193    /// \param cost The cost (length) values of the edges.
194    /// \param supply The supply values of the nodes (signed).
195    CycleCanceling( const Graph &graph,
196                    const CapacityMap &capacity,
197                    const CostMap &cost,
198                    const SupplyMap &supply ) :
199      _graph(graph), _lower(NULL), _capacity(capacity), _cost(cost),
200      _supply(supply), _flow(graph, 0),
201      _res_graph(graph, _capacity, _flow), _res_cost(_cost)
202    {
203      // Checking the sum of supply values
204      Supply sum = 0;
205      for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
206      _valid_supply = sum == 0;
207    }
208
209    /// \brief Simple constructor of the class (with lower bounds).
210    ///
211    /// Simple constructor of the class (with lower bounds).
212    ///
213    /// \param graph The directed graph the algorithm runs on.
214    /// \param lower The lower bounds of the edges.
215    /// \param capacity The capacities (upper bounds) of the edges.
216    /// \param cost The cost (length) values of the edges.
217    /// \param s The source node.
218    /// \param t The target node.
219    /// \param flow_value The required amount of flow from node \c s
220    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
221    CycleCanceling( const Graph &graph,
222                    const LowerMap &lower,
223                    const CapacityMap &capacity,
224                    const CostMap &cost,
225                    Node s, Node t,
226                    Supply flow_value ) :
227      _graph(graph), _lower(&lower), _capacity(graph), _cost(cost),
228      _supply(graph), _flow(graph, 0),
229      _res_graph(graph, _capacity, _flow), _res_cost(_cost)
230    {
231      // Removing non-zero lower bounds
232      _capacity = subMap(capacity, lower);
233      for (NodeIt n(_graph); n != INVALID; ++n) {
234        Supply sum = 0;
235        if (n == s) sum =  flow_value;
236        if (n == t) sum = -flow_value;
237        for (InEdgeIt e(_graph, n); e != INVALID; ++e)
238          sum += lower[e];
239        for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
240          sum -= lower[e];
241        _supply[n] = sum;
242      }
243      _valid_supply = true;
244    }
245
246    /// \brief Simple constructor of the class (without lower bounds).
247    ///
248    /// Simple constructor of the class (without lower bounds).
249    ///
250    /// \param graph The directed graph the algorithm runs on.
251    /// \param capacity The capacities (upper bounds) of the edges.
252    /// \param cost The cost (length) values of the edges.
253    /// \param s The source node.
254    /// \param t The target node.
255    /// \param flow_value The required amount of flow from node \c s
256    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
257    CycleCanceling( const Graph &graph,
258                    const CapacityMap &capacity,
259                    const CostMap &cost,
260                    Node s, Node t,
261                    Supply flow_value ) :
262      _graph(graph), _lower(NULL), _capacity(capacity), _cost(cost),
263      _supply(graph, 0), _flow(graph, 0),
264      _res_graph(graph, _capacity, _flow), _res_cost(_cost)
265    {
266      _supply[s] =  flow_value;
267      _supply[t] = -flow_value;
268      _valid_supply = true;
269    }
270
271    /// \brief Runs the algorithm.
272    ///
273    /// Runs the algorithm.
274    ///
275    /// \param min_mean_cc Set this parameter to \c true to run the
276    /// "Minimum Mean Cycle-Canceling" algorithm, which is strongly
277    /// polynomial, but rather slower in practice.
278    ///
279    /// \return \c true if a feasible flow can be found.
280    bool run(bool min_mean_cc = false) {
281      return init() && start(min_mean_cc);
282    }
283
284    /// \brief Returns a const reference to the edge map storing the
285    /// found flow.
286    ///
287    /// Returns a const reference to the edge map storing the found flow.
288    ///
289    /// \pre \ref run() must be called before using this function.
290    const FlowMap& flowMap() const {
291      return _flow;
292    }
293
294    /// \brief Returns the total cost of the found flow.
295    ///
296    /// Returns the total cost of the found flow. The complexity of the
297    /// function is \f$ O(e) \f$.
298    ///
299    /// \pre \ref run() must be called before using this function.
300    Cost totalCost() const {
301      Cost c = 0;
302      for (EdgeIt e(_graph); e != INVALID; ++e)
303        c += _flow[e] * _cost[e];
304      return c;
305    }
306
307  private:
308
309    /// Initializes the algorithm.
310    bool init() {
311      if (!_valid_supply) return false;
312
313      // Finding a feasible flow using Circulation
314      Circulation< Graph, ConstMap<Edge, Capacity>, CapacityEdgeMap,
315                   SupplyMap >
316        circulation( _graph, constMap<Edge>((Capacity)0), _capacity,
317                     _supply );
318      return circulation.flowMap(_flow).run();
319    }
320
321    bool start(bool min_mean_cc) {
322      if (min_mean_cc)
323        return startMinMean();
324      else
325        return start();
326    }
327
328    /// \brief Executes the algorithm using \ref BellmanFord.
329    ///
330    /// Executes the algorithm using the \ref BellmanFord
331    /// "Bellman-Ford" algorithm for negative cycle detection with
332    /// successively larger limit for the number of iterations.
333    bool start() {
334      typename BellmanFord<ResGraph, ResidualCostMap>::PredMap pred(_res_graph);
335      typename ResGraph::template NodeMap<int> visited(_res_graph);
336      std::vector<ResEdge> cycle;
337      int node_num = countNodes(_graph);
338
339      int length_bound = BF_FIRST_LIMIT;
340      bool optimal = false;
341      while (!optimal) {
342        BellmanFord<ResGraph, ResidualCostMap> bf(_res_graph, _res_cost);
343        bf.predMap(pred);
344        bf.init(0);
345        int iter_num = 0;
346        bool cycle_found = false;
347        while (!cycle_found) {
348          int curr_iter_num = iter_num + length_bound <= node_num ?
349                              length_bound : node_num - iter_num;
350          iter_num += curr_iter_num;
351          int real_iter_num = curr_iter_num;
352          for (int i = 0; i < curr_iter_num; ++i) {
353            if (bf.processNextWeakRound()) {
354              real_iter_num = i;
355              break;
356            }
357          }
358          if (real_iter_num < curr_iter_num) {
359            optimal = true;
360            break;
361          } else {
362            // Searching for node disjoint negative cycles
363            for (ResNodeIt n(_res_graph); n != INVALID; ++n)
364              visited[n] = 0;
365            int id = 0;
366            for (ResNodeIt n(_res_graph); n != INVALID; ++n) {
367              if (visited[n] > 0) continue;
368              visited[n] = ++id;
369              ResNode u = pred[n] == INVALID ?
370                          INVALID : _res_graph.source(pred[n]);
371              while (u != INVALID && visited[u] == 0) {
372                visited[u] = id;
373                u = pred[u] == INVALID ?
374                    INVALID : _res_graph.source(pred[u]);
375              }
376              if (u != INVALID && visited[u] == id) {
377                // Finding the negative cycle
378                cycle_found = true;
379                cycle.clear();
380                ResEdge e = pred[u];
381                cycle.push_back(e);
382                Capacity d = _res_graph.rescap(e);
383                while (_res_graph.source(e) != u) {
384                  cycle.push_back(e = pred[_res_graph.source(e)]);
385                  if (_res_graph.rescap(e) < d)
386                    d = _res_graph.rescap(e);
387                }
388
389                // Augmenting along the cycle
390                for (int i = 0; i < int(cycle.size()); ++i)
391                  _res_graph.augment(cycle[i], d);
392              }
393            }
394          }
395
396          if (!cycle_found)
397            length_bound = int(length_bound * BF_ALPHA);
398        }
399      }
400
401      // Handling non-zero lower bounds
402      if (_lower) {
403        for (EdgeIt e(_graph); e != INVALID; ++e)
404          _flow[e] += (*_lower)[e];
405      }
406      return true;
407    }
408
409    /// \brief Executes the algorithm using \ref MinMeanCycle.
410    ///
411    /// Executes the algorithm using \ref MinMeanCycle for negative
412    /// cycle detection.
413    bool startMinMean() {
414      typedef Path<ResGraph> ResPath;
415      MinMeanCycle<ResGraph, ResidualCostMap> mmc(_res_graph, _res_cost);
416      ResPath cycle;
417
418      mmc.cyclePath(cycle).init();
419      if (mmc.findMinMean()) {
420        while (mmc.cycleLength() < 0) {
421          // Finding the cycle
422          mmc.findCycle();
423
424          // Finding the largest flow amount that can be augmented
425          // along the cycle
426          Capacity delta = 0;
427          for (typename ResPath::EdgeIt e(cycle); e != INVALID; ++e) {
428            if (delta == 0 || _res_graph.rescap(e) < delta)
429              delta = _res_graph.rescap(e);
430          }
431
432          // Augmenting along the cycle
433          for (typename ResPath::EdgeIt e(cycle); e != INVALID; ++e)
434            _res_graph.augment(e, delta);
435
436          // Finding the minimum cycle mean for the modified residual
437          // graph
438          mmc.reset();
439          if (!mmc.findMinMean()) break;
440        }
441      }
442
443      // Handling non-zero lower bounds
444      if (_lower) {
445        for (EdgeIt e(_graph); e != INVALID; ++e)
446          _flow[e] += (*_lower)[e];
447      }
448      return true;
449    }
450
451  }; //class CycleCanceling
452
453  ///@}
454
455} //namespace lemon
456
457#endif //LEMON_CYCLE_CANCELING_H
Note: See TracBrowser for help on using the repository browser.