COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/cycle_canceling.h @ 2457:8c791ee69a45

Last change on this file since 2457:8c791ee69a45 was 2457:8c791ee69a45, checked in by Balazs Dezso, 12 years ago

Improvments in min cost flow algorithms

  • improved cycle cancelling
File size: 14.9 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_CYCLE_CANCELING_H
20#define LEMON_CYCLE_CANCELING_H
21
22/// \ingroup min_cost_flow
23///
24/// \file
25/// \brief A cycle-canceling algorithm for finding a minimum cost flow.
26
27#include <vector>
28#include <lemon/circulation.h>
29#include <lemon/graph_adaptor.h>
30
31/// \brief The used cycle-canceling method.
32#define LIMITED_CYCLE_CANCELING
33//#define MIN_MEAN_CYCLE_CANCELING
34
35#ifdef LIMITED_CYCLE_CANCELING
36  #include <lemon/bellman_ford.h>
37  /// \brief The maximum number of iterations for the first execution
38  /// of the \ref lemon::BellmanFord "Bellman-Ford" algorithm.
39  /// It should be at least 2.
40  #define STARTING_LIMIT        2
41  /// \brief The iteration limit for the
42  /// \ref lemon::BellmanFord "Bellman-Ford" algorithm is multiplied by
43  /// <tt>ALPHA_MUL / ALPHA_DIV</tt> in every round.
44  /// <tt>ALPHA_MUL / ALPHA_DIV</tt> must be greater than 1.
45  #define ALPHA_MUL             3
46  /// \brief The iteration limit for the
47  /// \ref lemon::BellmanFord "Bellman-Ford" algorithm is multiplied by
48  /// <tt>ALPHA_MUL / ALPHA_DIV</tt> in every round.
49  /// <tt>ALPHA_MUL / ALPHA_DIV</tt> must be greater than 1.
50  #define ALPHA_DIV             2
51
52//#define _ONLY_ONE_CYCLE_
53//#define _NO_BACK_STEP_
54//#define _DEBUG_ITER_
55#endif
56
57#ifdef MIN_MEAN_CYCLE_CANCELING
58  #include <lemon/min_mean_cycle.h>
59  #include <lemon/path.h>
60#endif
61
62namespace lemon {
63
64  /// \addtogroup min_cost_flow
65  /// @{
66
67  /// \brief Implementation of a cycle-canceling algorithm for finding
68  /// a minimum cost flow.
69  ///
70  /// \ref lemon::CycleCanceling "CycleCanceling" implements a
71  /// cycle-canceling algorithm for finding a minimum cost flow.
72  ///
73  /// \param Graph The directed graph type the algorithm runs on.
74  /// \param LowerMap The type of the lower bound map.
75  /// \param CapacityMap The type of the capacity (upper bound) map.
76  /// \param CostMap The type of the cost (length) map.
77  /// \param SupplyMap The type of the supply map.
78  ///
79  /// \warning
80  /// - Edge capacities and costs should be nonnegative integers.
81  ///   However \c CostMap::Value should be signed type.
82  /// - Supply values should be integers.
83  /// - \c LowerMap::Value must be convertible to
84  ///   \c CapacityMap::Value and \c CapacityMap::Value must be
85  ///   convertible to \c SupplyMap::Value.
86  ///
87  /// \author Peter Kovacs
88
89template < typename Graph,
90           typename LowerMap = typename Graph::template EdgeMap<int>,
91           typename CapacityMap = LowerMap,
92           typename CostMap = typename Graph::template EdgeMap<int>,
93           typename SupplyMap = typename Graph::template NodeMap
94                                <typename CapacityMap::Value> >
95  class CycleCanceling
96  {
97    typedef typename Graph::Node Node;
98    typedef typename Graph::NodeIt NodeIt;
99    typedef typename Graph::Edge Edge;
100    typedef typename Graph::EdgeIt EdgeIt;
101    typedef typename Graph::InEdgeIt InEdgeIt;
102    typedef typename Graph::OutEdgeIt OutEdgeIt;
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    typedef typename Graph::template EdgeMap<Capacity> CapacityRefMap;
109    typedef typename Graph::template NodeMap<Supply> SupplyRefMap;
110
111    typedef ResGraphAdaptor< const Graph, Capacity,
112                             CapacityRefMap, CapacityRefMap > ResGraph;
113    typedef typename ResGraph::Node ResNode;
114    typedef typename ResGraph::NodeIt ResNodeIt;
115    typedef typename ResGraph::Edge ResEdge;
116    typedef typename ResGraph::EdgeIt ResEdgeIt;
117
118  public:
119
120    /// \brief The type of the flow map.
121    typedef CapacityRefMap FlowMap;
122
123  protected:
124
125    /// \brief Map adaptor class for handling residual edge costs.
126    class ResCostMap : public MapBase<ResEdge, Cost>
127    {
128    private:
129
130      const CostMap &cost_map;
131
132    public:
133
134      typedef typename MapBase<ResEdge, Cost>::Value Value;
135      typedef typename MapBase<ResEdge, Cost>::Key Key;
136
137      ResCostMap(const CostMap &_cost) : cost_map(_cost) {}
138
139      Value operator[](const Key &e) const {
140        return ResGraph::forward(e) ? cost_map[e] : -cost_map[e];
141      }
142
143    }; //class ResCostMap
144
145  protected:
146
147    /// \brief The directed graph the algorithm runs on.
148    const Graph &graph;
149    /// \brief The original lower bound map.
150    const LowerMap *lower;
151    /// \brief The modified capacity map.
152    CapacityRefMap capacity;
153    /// \brief The cost map.
154    const CostMap &cost;
155    /// \brief The modified supply map.
156    SupplyRefMap supply;
157    /// \brief The sum of supply values equals zero.
158    bool valid_supply;
159
160    /// \brief The current flow.
161    FlowMap flow;
162    /// \brief The residual graph.
163    ResGraph res_graph;
164    /// \brief The residual cost map.
165    ResCostMap res_cost;
166
167  public :
168
169    /// \brief General constructor of the class (with lower bounds).
170    ///
171    /// General constructor of the class (with lower bounds).
172    ///
173    /// \param _graph The directed graph the algorithm runs on.
174    /// \param _lower The lower bounds of the edges.
175    /// \param _capacity The capacities (upper bounds) of the edges.
176    /// \param _cost The cost (length) values of the edges.
177    /// \param _supply The supply values of the nodes (signed).
178    CycleCanceling( const Graph &_graph,
179                    const LowerMap &_lower,
180                    const CapacityMap &_capacity,
181                    const CostMap &_cost,
182                    const SupplyMap &_supply ) :
183      graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
184      supply(_graph), flow(_graph, 0),
185      res_graph(_graph, capacity, flow), res_cost(cost)
186    {
187      // Removing nonzero lower bounds
188      capacity = subMap(_capacity, _lower);
189      Supply sum = 0;
190      for (NodeIt n(graph); n != INVALID; ++n) {
191        Supply s = _supply[n];
192        for (InEdgeIt e(graph, n); e != INVALID; ++e)
193          s += _lower[e];
194        for (OutEdgeIt e(graph, n); e != INVALID; ++e)
195          s -= _lower[e];
196        sum += (supply[n] = s);
197      }
198      valid_supply = sum == 0;
199    }
200
201    /// \brief General constructor of the class (without lower bounds).
202    ///
203    /// General constructor of the class (without lower bounds).
204    ///
205    /// \param _graph The directed graph the algorithm runs on.
206    /// \param _capacity The capacities (upper bounds) of the edges.
207    /// \param _cost The cost (length) values of the edges.
208    /// \param _supply The supply values of the nodes (signed).
209    CycleCanceling( const Graph &_graph,
210                    const CapacityMap &_capacity,
211                    const CostMap &_cost,
212                    const SupplyMap &_supply ) :
213      graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
214      supply(_supply), flow(_graph, 0),
215      res_graph(_graph, capacity, flow), res_cost(cost)
216    {
217      // Checking the sum of supply values
218      Supply sum = 0;
219      for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n];
220      valid_supply = sum == 0;
221    }
222
223
224    /// \brief Simple constructor of the class (with lower bounds).
225    ///
226    /// Simple constructor of the class (with lower bounds).
227    ///
228    /// \param _graph The directed graph the algorithm runs on.
229    /// \param _lower The lower bounds of the edges.
230    /// \param _capacity The capacities (upper bounds) of the edges.
231    /// \param _cost The cost (length) values of the edges.
232    /// \param _s The source node.
233    /// \param _t The target node.
234    /// \param _flow_value The required amount of flow from node \c _s
235    /// to node \c _t (i.e. the supply of \c _s and the demand of
236    /// \c _t).
237    CycleCanceling( const Graph &_graph,
238                    const LowerMap &_lower,
239                    const CapacityMap &_capacity,
240                    const CostMap &_cost,
241                    Node _s, Node _t,
242                    Supply _flow_value ) :
243      graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
244      supply(_graph), flow(_graph, 0),
245      res_graph(_graph, capacity, flow), res_cost(cost)
246    {
247      // Removing nonzero lower bounds
248      capacity = subMap(_capacity, _lower);
249      for (NodeIt n(graph); n != INVALID; ++n) {
250        Supply s = 0;
251        if (n == _s) s =  _flow_value;
252        if (n == _t) s = -_flow_value;
253        for (InEdgeIt e(graph, n); e != INVALID; ++e)
254          s += _lower[e];
255        for (OutEdgeIt e(graph, n); e != INVALID; ++e)
256          s -= _lower[e];
257        supply[n] = s;
258      }
259      valid_supply = true;
260    }
261
262    /// \brief Simple constructor of the class (without lower bounds).
263    ///
264    /// Simple constructor of the class (without lower bounds).
265    ///
266    /// \param _graph The directed graph the algorithm runs on.
267    /// \param _capacity The capacities (upper bounds) of the edges.
268    /// \param _cost The cost (length) values of the edges.
269    /// \param _s The source node.
270    /// \param _t The target node.
271    /// \param _flow_value The required amount of flow from node \c _s
272    /// to node \c _t (i.e. the supply of \c _s and the demand of
273    /// \c _t).
274    CycleCanceling( const Graph &_graph,
275                    const CapacityMap &_capacity,
276                    const CostMap &_cost,
277                    Node _s, Node _t,
278                    Supply _flow_value ) :
279      graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
280      supply(_graph, 0), flow(_graph, 0),
281      res_graph(_graph, capacity, flow), res_cost(cost)
282    {
283      supply[_s] =  _flow_value;
284      supply[_t] = -_flow_value;
285      valid_supply = true;
286    }
287
288    /// \brief Returns a const reference to the flow map.
289    ///
290    /// Returns a const reference to the flow map.
291    ///
292    /// \pre \ref run() must be called before using this function.
293    const FlowMap& flowMap() const {
294      return flow;
295    }
296
297    /// \brief Returns the total cost of the found flow.
298    ///
299    /// Returns the total cost of the found flow. The complexity of the
300    /// function is \f$ O(e) \f$.
301    ///
302    /// \pre \ref run() must be called before using this function.
303    Cost totalCost() const {
304      Cost c = 0;
305      for (EdgeIt e(graph); e != INVALID; ++e)
306        c += flow[e] * cost[e];
307      return c;
308    }
309
310    /// \brief Runs the algorithm.
311    ///
312    /// Runs the algorithm.
313    ///
314    /// \return \c true if a feasible flow can be found.
315    bool run() {
316      return init() && start();
317    }
318
319  protected:
320
321    /// \brief Initializes the algorithm.
322    bool init() {
323      // Checking the sum of supply values
324      Supply sum = 0;
325      for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n];
326      if (sum != 0) return false;
327
328      // Finding a feasible flow
329      Circulation< Graph, Capacity, FlowMap, ConstMap<Edge, Capacity>,
330                   CapacityRefMap, SupplyMap >
331        circulation( graph, constMap<Edge>((Capacity)0),
332                     capacity, supply, flow );
333      return circulation.run() == -1;
334    }
335
336#ifdef LIMITED_CYCLE_CANCELING
337    /// \brief Executes a cycle-canceling algorithm using
338    /// \ref lemon::BellmanFord "Bellman-Ford" algorithm with limited
339    /// iteration count.
340    bool start() {
341      typename BellmanFord<ResGraph, ResCostMap>::PredMap pred(res_graph);
342      typename ResGraph::template NodeMap<int> visited(res_graph);
343      std::vector<ResEdge> cycle;
344      int node_num = countNodes(graph);
345
346#ifdef _DEBUG_ITER_
347      int cycle_num = 0;
348#endif
349      int length_bound = STARTING_LIMIT;
350      bool optimal = false;
351      while (!optimal) {
352        BellmanFord<ResGraph, ResCostMap> bf(res_graph, res_cost);
353        bf.predMap(pred);
354        bf.init(0);
355        int iter_num = 0;
356        bool cycle_found = false;
357        while (!cycle_found) {
358#ifdef _NO_BACK_STEP_
359          int curr_iter_num = length_bound <= node_num ?
360                              length_bound - iter_num : node_num - iter_num;
361#else
362          int curr_iter_num = iter_num + length_bound <= node_num ?
363                              length_bound : node_num - iter_num;
364#endif
365          iter_num += curr_iter_num;
366          int real_iter_num = curr_iter_num;
367          for (int i = 0; i < curr_iter_num; ++i) {
368            if (bf.processNextWeakRound()) {
369              real_iter_num = i;
370              break;
371            }
372          }
373          if (real_iter_num < curr_iter_num) {
374            optimal = true;
375            break;
376          } else {
377            // Searching for node disjoint negative cycles
378            for (ResNodeIt n(res_graph); n != INVALID; ++n)
379              visited[n] = 0;
380            int id = 0;
381            for (ResNodeIt n(res_graph); n != INVALID; ++n) {
382              if (visited[n] > 0) continue;
383              visited[n] = ++id;
384              ResNode u = pred[n] == INVALID ?
385                          INVALID : res_graph.source(pred[n]);
386              while (u != INVALID && visited[u] == 0) {
387                visited[u] = id;
388                u = pred[u] == INVALID ?
389                    INVALID : res_graph.source(pred[u]);
390              }
391              if (u != INVALID && visited[u] == id) {
392                // Finding the negative cycle
393                cycle_found = true;
394                cycle.clear();
395                ResEdge e = pred[u];
396                cycle.push_back(e);
397                Capacity d = res_graph.rescap(e);
398                while (res_graph.source(e) != u) {
399                  cycle.push_back(e = pred[res_graph.source(e)]);
400                  if (res_graph.rescap(e) < d)
401                    d = res_graph.rescap(e);
402                }
403#ifdef _DEBUG_ITER_
404                ++cycle_num;
405#endif
406                // Augmenting along the cycle
407                for (int i = 0; i < cycle.size(); ++i)
408                  res_graph.augment(cycle[i], d);
409#ifdef _ONLY_ONE_CYCLE_
410                break;
411#endif
412              }
413            }
414          }
415
416          if (!cycle_found)
417            length_bound = length_bound * ALPHA_MUL / ALPHA_DIV;
418        }
419      }
420
421#ifdef _DEBUG_ITER_
422      std::cout << "Limited cycle-canceling algorithm finished. "
423                << "Found " << cycle_num << " negative cycles."
424                << std::endl;
425#endif
426
427      // Handling nonzero lower bounds
428      if (lower) {
429        for (EdgeIt e(graph); e != INVALID; ++e)
430          flow[e] += (*lower)[e];
431      }
432      return true;
433    }
434#endif
435
436#ifdef MIN_MEAN_CYCLE_CANCELING
437    /// \brief Executes the minimum mean cycle-canceling algorithm
438    /// using \ref lemon::MinMeanCycle "MinMeanCycle" class.
439    bool start() {
440      typedef Path<ResGraph> ResPath;
441      MinMeanCycle<ResGraph, ResCostMap> mmc(res_graph, res_cost);
442      ResPath cycle;
443
444#ifdef _DEBUG_ITER_
445      int cycle_num = 0;
446#endif
447      mmc.cyclePath(cycle).init();
448      if (mmc.findMinMean()) {
449        while (mmc.cycleLength() < 0) {
450#ifdef _DEBUG_ITER_
451          ++iter;
452#endif
453          // Finding the cycle
454          mmc.findCycle();
455
456          // Finding the largest flow amount that can be augmented
457          // along the cycle
458          Capacity delta = 0;
459          for (typename ResPath::EdgeIt e(cycle); e != INVALID; ++e) {
460            if (delta == 0 || res_graph.rescap(e) < delta)
461              delta = res_graph.rescap(e);
462          }
463
464          // Augmenting along the cycle
465          for (typename ResPath::EdgeIt e(cycle); e != INVALID; ++e)
466            res_graph.augment(e, delta);
467
468          // Finding the minimum cycle mean for the modified residual
469          // graph
470          mmc.reset();
471          if (!mmc.findMinMean()) break;
472        }
473      }
474
475#ifdef _DEBUG_ITER_
476      std::cout << "Minimum mean cycle-canceling algorithm finished. "
477                << "Found " << cycle_num << " negative cycles."
478                << std::endl;
479#endif
480
481      // Handling nonzero lower bounds
482      if (lower) {
483        for (EdgeIt e(graph); e != INVALID; ++e)
484          flow[e] += (*lower)[e];
485      }
486      return true;
487    }
488#endif
489
490  }; //class CycleCanceling
491
492  ///@}
493
494} //namespace lemon
495
496#endif //LEMON_CYCLE_CANCELING_H
Note: See TracBrowser for help on using the repository browser.