COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/capacity_scaling.h @ 2574:7058c9690e7d

Last change on this file since 2574:7058c9690e7d was 2574:7058c9690e7d, checked in by Peter Kovacs, 14 years ago

Improvements in CapacityScaling?.

Main changes:

  • Use -potenital[] instead of potential[] to conform to the usual

terminology.

  • 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: 20.1 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_CAPACITY_SCALING_H
20#define LEMON_CAPACITY_SCALING_H
21
22/// \ingroup min_cost_flow
23///
24/// \file
25/// \brief Capacity scaling algorithm for finding a minimum cost flow.
26
27#include <vector>
28
29#include <lemon/graph_adaptor.h>
30#include <lemon/bin_heap.h>
31
32namespace lemon {
33
34  /// \addtogroup min_cost_flow
35  /// @{
36
37  /// \brief Implementation of the capacity scaling algorithm for
38  /// finding a minimum cost flow.
39  ///
40  /// \ref CapacityScaling implements the capacity scaling version
41  /// of the successive shortest path algorithm for finding a minimum
42  /// cost flow.
43  ///
44  /// \tparam Graph The directed graph type the algorithm runs on.
45  /// \tparam LowerMap The type of the lower bound map.
46  /// \tparam CapacityMap The type of the capacity (upper bound) map.
47  /// \tparam CostMap The type of the cost (length) map.
48  /// \tparam SupplyMap The type of the supply map.
49  ///
50  /// \warning
51  /// - Edge capacities and costs should be \e non-negative \e integers.
52  /// - Supply values should be \e signed \e integers.
53  /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.
54  /// - \c CapacityMap::Value and \c SupplyMap::Value must be
55  ///   convertible to each other.
56  /// - All value types must be convertible to \c CostMap::Value, which
57  ///   must be signed type.
58  ///
59  /// \author Peter Kovacs
60
61  template < typename Graph,
62             typename LowerMap = typename Graph::template EdgeMap<int>,
63             typename CapacityMap = typename Graph::template EdgeMap<int>,
64             typename CostMap = typename Graph::template EdgeMap<int>,
65             typename SupplyMap = typename Graph::template NodeMap<int> >
66  class CapacityScaling
67  {
68    GRAPH_TYPEDEFS(typename Graph);
69
70    typedef typename CapacityMap::Value Capacity;
71    typedef typename CostMap::Value Cost;
72    typedef typename SupplyMap::Value Supply;
73    typedef typename Graph::template EdgeMap<Capacity> CapacityEdgeMap;
74    typedef typename Graph::template NodeMap<Supply> SupplyNodeMap;
75    typedef typename Graph::template NodeMap<Edge> PredMap;
76
77  public:
78
79    /// The type of the flow map.
80    typedef typename Graph::template EdgeMap<Capacity> FlowMap;
81    /// The type of the potential map.
82    typedef typename Graph::template NodeMap<Cost> PotentialMap;
83
84  private:
85
86    /// \brief Special implementation of the \ref Dijkstra algorithm
87    /// for finding shortest paths in the residual network.
88    ///
89    /// \ref ResidualDijkstra is a special implementation of the
90    /// \ref Dijkstra algorithm for finding shortest paths in the
91    /// residual network of the graph with respect to the reduced edge
92    /// costs and modifying the node potentials according to the
93    /// distance of the nodes.
94    class ResidualDijkstra
95    {
96      typedef typename Graph::template NodeMap<Cost> CostNodeMap;
97      typedef typename Graph::template NodeMap<Edge> PredMap;
98
99      typedef typename Graph::template NodeMap<int> HeapCrossRef;
100      typedef BinHeap<Cost, HeapCrossRef> Heap;
101
102    private:
103
104      // The directed graph the algorithm runs on
105      const Graph &_graph;
106
107      // The main maps
108      const FlowMap &_flow;
109      const CapacityEdgeMap &_res_cap;
110      const CostMap &_cost;
111      const SupplyNodeMap &_excess;
112      PotentialMap &_potential;
113
114      // The distance map
115      CostNodeMap _dist;
116      // The pred edge map
117      PredMap &_pred;
118      // The processed (i.e. permanently labeled) nodes
119      std::vector<Node> _proc_nodes;
120
121    public:
122
123      /// The constructor of the class.
124      ResidualDijkstra( const Graph &graph,
125                        const FlowMap &flow,
126                        const CapacityEdgeMap &res_cap,
127                        const CostMap &cost,
128                        const SupplyMap &excess,
129                        PotentialMap &potential,
130                        PredMap &pred ) :
131        _graph(graph), _flow(flow), _res_cap(res_cap), _cost(cost),
132        _excess(excess), _potential(potential), _dist(graph),
133        _pred(pred)
134      {}
135
136      /// Runs the algorithm from the given source node.
137      Node run(Node s, Capacity delta) {
138        HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
139        Heap heap(heap_cross_ref);
140        heap.push(s, 0);
141        _pred[s] = INVALID;
142        _proc_nodes.clear();
143
144        // Processing nodes
145        while (!heap.empty() && _excess[heap.top()] > -delta) {
146          Node u = heap.top(), v;
147          Cost d = heap.prio() + _potential[u], nd;
148          _dist[u] = heap.prio();
149          heap.pop();
150          _proc_nodes.push_back(u);
151
152          // Traversing outgoing edges
153          for (OutEdgeIt e(_graph, u); e != INVALID; ++e) {
154            if (_res_cap[e] >= delta) {
155              v = _graph.target(e);
156              switch(heap.state(v)) {
157              case Heap::PRE_HEAP:
158                heap.push(v, d + _cost[e] - _potential[v]);
159                _pred[v] = e;
160                break;
161              case Heap::IN_HEAP:
162                nd = d + _cost[e] - _potential[v];
163                if (nd < heap[v]) {
164                  heap.decrease(v, nd);
165                  _pred[v] = e;
166                }
167                break;
168              case Heap::POST_HEAP:
169                break;
170              }
171            }
172          }
173
174          // Traversing incoming edges
175          for (InEdgeIt e(_graph, u); e != INVALID; ++e) {
176            if (_flow[e] >= delta) {
177              v = _graph.source(e);
178              switch(heap.state(v)) {
179              case Heap::PRE_HEAP:
180                heap.push(v, d - _cost[e] - _potential[v]);
181                _pred[v] = e;
182                break;
183              case Heap::IN_HEAP:
184                nd = d - _cost[e] - _potential[v];
185                if (nd < heap[v]) {
186                  heap.decrease(v, nd);
187                  _pred[v] = e;
188                }
189                break;
190              case Heap::POST_HEAP:
191                break;
192              }
193            }
194          }
195        }
196        if (heap.empty()) return INVALID;
197
198        // Updating potentials of processed nodes
199        Node t = heap.top();
200        Cost t_dist = heap.prio();
201        for (int i = 0; i < int(_proc_nodes.size()); ++i)
202          _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
203
204        return t;
205      }
206
207    }; //class ResidualDijkstra
208
209  private:
210
211    // The directed graph the algorithm runs on
212    const Graph &_graph;
213    // The original lower bound map
214    const LowerMap *_lower;
215    // The modified capacity map
216    CapacityEdgeMap _capacity;
217    // The original cost map
218    const CostMap &_cost;
219    // The modified supply map
220    SupplyNodeMap _supply;
221    bool _valid_supply;
222
223    // Edge map of the current flow
224    FlowMap _flow;
225    // Node map of the current potentials
226    PotentialMap _potential;
227
228    // The residual capacity map
229    CapacityEdgeMap _res_cap;
230    // The excess map
231    SupplyNodeMap _excess;
232    // The excess nodes (i.e. nodes with positive excess)
233    std::vector<Node> _excess_nodes;
234    // The deficit nodes (i.e. nodes with negative excess)
235    std::vector<Node> _deficit_nodes;
236
237    // The delta parameter used for capacity scaling
238    Capacity _delta;
239    // The maximum number of phases
240    int _phase_num;
241
242    // The pred edge map
243    PredMap _pred;
244    // Implementation of the Dijkstra algorithm for finding augmenting
245    // shortest paths in the residual network
246    ResidualDijkstra _dijkstra;
247
248  public :
249
250    /// \brief General constructor of the class (with lower bounds).
251    ///
252    /// General constructor of the class (with lower bounds).
253    ///
254    /// \param graph The directed graph the algorithm runs on.
255    /// \param lower The lower bounds of the edges.
256    /// \param capacity The capacities (upper bounds) of the edges.
257    /// \param cost The cost (length) values of the edges.
258    /// \param supply The supply values of the nodes (signed).
259    CapacityScaling( const Graph &graph,
260                     const LowerMap &lower,
261                     const CapacityMap &capacity,
262                     const CostMap &cost,
263                     const SupplyMap &supply ) :
264      _graph(graph), _lower(&lower), _capacity(graph), _cost(cost),
265      _supply(graph), _flow(graph, 0), _potential(graph, 0),
266      _res_cap(graph), _excess(graph), _pred(graph),
267      _dijkstra(_graph, _flow, _res_cap, _cost, _excess, _potential, _pred)
268    {
269      // Removing non-zero lower bounds
270      _capacity = subMap(capacity, lower);
271      _res_cap = _capacity;
272      Supply sum = 0;
273      for (NodeIt n(_graph); n != INVALID; ++n) {
274        Supply s = supply[n];
275        for (InEdgeIt e(_graph, n); e != INVALID; ++e)
276          s += lower[e];
277        for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
278          s -= lower[e];
279        _supply[n] = s;
280        sum += s;
281      }
282      _valid_supply = sum == 0;
283    }
284
285    /// \brief General constructor of the class (without lower bounds).
286    ///
287    /// General constructor of the class (without lower bounds).
288    ///
289    /// \param graph The directed graph the algorithm runs on.
290    /// \param capacity The capacities (upper bounds) of the edges.
291    /// \param cost The cost (length) values of the edges.
292    /// \param supply The supply values of the nodes (signed).
293    CapacityScaling( const Graph &graph,
294                     const CapacityMap &capacity,
295                     const CostMap &cost,
296                     const SupplyMap &supply ) :
297      _graph(graph), _lower(NULL), _capacity(capacity), _cost(cost),
298      _supply(supply), _flow(graph, 0), _potential(graph, 0),
299      _res_cap(capacity), _excess(graph), _pred(graph),
300      _dijkstra(_graph, _flow, _res_cap, _cost, _excess, _potential, _pred)
301    {
302      // Checking the sum of supply values
303      Supply sum = 0;
304      for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
305      _valid_supply = sum == 0;
306    }
307
308    /// \brief Simple constructor of the class (with lower bounds).
309    ///
310    /// Simple constructor of the class (with lower bounds).
311    ///
312    /// \param graph The directed graph the algorithm runs on.
313    /// \param lower The lower bounds of the edges.
314    /// \param capacity The capacities (upper bounds) of the edges.
315    /// \param cost The cost (length) values of the edges.
316    /// \param s The source node.
317    /// \param t The target node.
318    /// \param flow_value The required amount of flow from node \c s
319    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
320    CapacityScaling( const Graph &graph,
321                     const LowerMap &lower,
322                     const CapacityMap &capacity,
323                     const CostMap &cost,
324                     Node s, Node t,
325                     Supply flow_value ) :
326      _graph(graph), _lower(&lower), _capacity(graph), _cost(cost),
327      _supply(graph), _flow(graph, 0), _potential(graph, 0),
328      _res_cap(graph), _excess(graph), _pred(graph),
329      _dijkstra(_graph, _flow, _res_cap, _cost, _excess, _potential, _pred)
330    {
331      // Removing non-zero lower bounds
332      _capacity = subMap(capacity, lower);
333      _res_cap = _capacity;
334      for (NodeIt n(_graph); n != INVALID; ++n) {
335        Supply sum = 0;
336        if (n == s) sum =  flow_value;
337        if (n == t) sum = -flow_value;
338        for (InEdgeIt e(_graph, n); e != INVALID; ++e)
339          sum += lower[e];
340        for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
341          sum -= lower[e];
342        _supply[n] = sum;
343      }
344      _valid_supply = true;
345    }
346
347    /// \brief Simple constructor of the class (without lower bounds).
348    ///
349    /// Simple constructor of the class (without lower bounds).
350    ///
351    /// \param graph The directed graph the algorithm runs on.
352    /// \param capacity The capacities (upper bounds) of the edges.
353    /// \param cost The cost (length) values of the edges.
354    /// \param s The source node.
355    /// \param t The target node.
356    /// \param flow_value The required amount of flow from node \c s
357    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
358    CapacityScaling( const Graph &graph,
359                     const CapacityMap &capacity,
360                     const CostMap &cost,
361                     Node s, Node t,
362                     Supply flow_value ) :
363      _graph(graph), _lower(NULL), _capacity(capacity), _cost(cost),
364      _supply(graph, 0), _flow(graph, 0), _potential(graph, 0),
365      _res_cap(capacity), _excess(graph), _pred(graph),
366      _dijkstra(_graph, _flow, _res_cap, _cost, _excess, _potential, _pred)
367    {
368      _supply[s] =  flow_value;
369      _supply[t] = -flow_value;
370      _valid_supply = true;
371    }
372
373    /// \brief Runs the algorithm.
374    ///
375    /// Runs the algorithm.
376    ///
377    /// \param scaling Enable or disable capacity scaling.
378    /// If the maximum edge capacity and/or the amount of total supply
379    /// is rather small, the algorithm could be slightly faster without
380    /// scaling.
381    ///
382    /// \return \c true if a feasible flow can be found.
383    bool run(bool scaling = true) {
384      return init(scaling) && start();
385    }
386
387    /// \brief Returns a const reference to the edge map storing the
388    /// found flow.
389    ///
390    /// Returns a const reference to the edge map storing the found flow.
391    ///
392    /// \pre \ref run() must be called before using this function.
393    const FlowMap& flowMap() const {
394      return _flow;
395    }
396
397    /// \brief Returns a const reference to the node map storing the
398    /// found potentials (the dual solution).
399    ///
400    /// Returns a const reference to the node map storing the found
401    /// potentials (the dual solution).
402    ///
403    /// \pre \ref run() must be called before using this function.
404    const PotentialMap& potentialMap() const {
405      return _potential;
406    }
407
408    /// \brief Returns the total cost of the found flow.
409    ///
410    /// Returns the total cost of the found flow. The complexity of the
411    /// function is \f$ O(e) \f$.
412    ///
413    /// \pre \ref run() must be called before using this function.
414    Cost totalCost() const {
415      Cost c = 0;
416      for (EdgeIt e(_graph); e != INVALID; ++e)
417        c += _flow[e] * _cost[e];
418      return c;
419    }
420
421  private:
422
423    /// Initializes the algorithm.
424    bool init(bool scaling) {
425      if (!_valid_supply) return false;
426      _excess = _supply;
427
428      // Initilaizing delta value
429      if (scaling) {
430        // With scaling
431        Supply max_sup = 0, max_dem = 0;
432        for (NodeIt n(_graph); n != INVALID; ++n) {
433          if ( _supply[n] > max_sup) max_sup =  _supply[n];
434          if (-_supply[n] > max_dem) max_dem = -_supply[n];
435        }
436        if (max_dem < max_sup) max_sup = max_dem;
437        _phase_num = 0;
438        for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2)
439          ++_phase_num;
440      } else {
441        // Without scaling
442        _delta = 1;
443      }
444      return true;
445    }
446
447    bool start() {
448      if (_delta > 1)
449        return startWithScaling();
450      else
451        return startWithoutScaling();
452    }
453
454    /// Executes the capacity scaling algorithm.
455    bool startWithScaling() {
456      // Processing capacity scaling phases
457      Node s, t;
458      int phase_cnt = 0;
459      int factor = 4;
460      while (true) {
461        // Saturating all edges not satisfying the optimality condition
462        for (EdgeIt e(_graph); e != INVALID; ++e) {
463          Node u = _graph.source(e), v = _graph.target(e);
464          Cost c = _cost[e] + _potential[u] - _potential[v];
465          if (c < 0 && _res_cap[e] >= _delta) {
466            _excess[u] -= _res_cap[e];
467            _excess[v] += _res_cap[e];
468            _flow[e] = _capacity[e];
469            _res_cap[e] = 0;
470          }
471          else if (c > 0 && _flow[e] >= _delta) {
472            _excess[u] += _flow[e];
473            _excess[v] -= _flow[e];
474            _flow[e] = 0;
475            _res_cap[e] = _capacity[e];
476          }
477        }
478
479        // Finding excess nodes and deficit nodes
480        _excess_nodes.clear();
481        _deficit_nodes.clear();
482        for (NodeIt n(_graph); n != INVALID; ++n) {
483          if (_excess[n] >=  _delta) _excess_nodes.push_back(n);
484          if (_excess[n] <= -_delta) _deficit_nodes.push_back(n);
485        }
486        int next_node = 0;
487
488        // Finding augmenting shortest paths
489        while (next_node < int(_excess_nodes.size())) {
490          // Checking deficit nodes
491          if (_delta > 1) {
492            bool delta_deficit = false;
493            for (int i = 0; i < int(_deficit_nodes.size()); ++i) {
494              if (_excess[_deficit_nodes[i]] <= -_delta) {
495                delta_deficit = true;
496                break;
497              }
498            }
499            if (!delta_deficit) break;
500          }
501
502          // Running Dijkstra
503          s = _excess_nodes[next_node];
504          if ((t = _dijkstra.run(s, _delta)) == INVALID) {
505            if (_delta > 1) {
506              ++next_node;
507              continue;
508            }
509            return false;
510          }
511
512          // Augmenting along a shortest path from s to t.
513          Capacity d = _excess[s] < -_excess[t] ? _excess[s] : -_excess[t];
514          Node u = t;
515          Edge e;
516          if (d > _delta) {
517            while ((e = _pred[u]) != INVALID) {
518              Capacity rc;
519              if (u == _graph.target(e)) {
520                rc = _res_cap[e];
521                u = _graph.source(e);
522              } else {
523                rc = _flow[e];
524                u = _graph.target(e);
525              }
526              if (rc < d) d = rc;
527            }
528          }
529          u = t;
530          while ((e = _pred[u]) != INVALID) {
531            if (u == _graph.target(e)) {
532              _flow[e] += d;
533              _res_cap[e] -= d;
534              u = _graph.source(e);
535            } else {
536              _flow[e] -= d;
537              _res_cap[e] += d;
538              u = _graph.target(e);
539            }
540          }
541          _excess[s] -= d;
542          _excess[t] += d;
543
544          if (_excess[s] < _delta) ++next_node;
545        }
546
547        if (_delta == 1) break;
548        if (++phase_cnt > _phase_num / 4) factor = 2;
549        _delta = _delta <= factor ? 1 : _delta / factor;
550      }
551
552      // Handling non-zero lower bounds
553      if (_lower) {
554        for (EdgeIt e(_graph); e != INVALID; ++e)
555          _flow[e] += (*_lower)[e];
556      }
557      return true;
558    }
559
560    /// Executes the successive shortest path algorithm.
561    bool startWithoutScaling() {
562      // Finding excess nodes
563      for (NodeIt n(_graph); n != INVALID; ++n)
564        if (_excess[n] > 0) _excess_nodes.push_back(n);
565      if (_excess_nodes.size() == 0) return true;
566      int next_node = 0;
567
568      // Finding shortest paths
569      Node s, t;
570      while ( _excess[_excess_nodes[next_node]] > 0 ||
571              ++next_node < int(_excess_nodes.size()) )
572      {
573        // Running Dijkstra
574        s = _excess_nodes[next_node];
575        if ((t = _dijkstra.run(s, 1)) == INVALID)
576          return false;
577
578        // Augmenting along a shortest path from s to t
579        Capacity d = _excess[s] < -_excess[t] ? _excess[s] : -_excess[t];
580        Node u = t;
581        Edge e;
582        while ((e = _pred[u]) != INVALID) {
583          Capacity rc;
584          if (u == _graph.target(e)) {
585            rc = _res_cap[e];
586            u = _graph.source(e);
587          } else {
588            rc = _flow[e];
589            u = _graph.target(e);
590          }
591          if (rc < d) d = rc;
592        }
593        u = t;
594        while ((e = _pred[u]) != INVALID) {
595          if (u == _graph.target(e)) {
596            _flow[e] += d;
597            _res_cap[e] -= d;
598            u = _graph.source(e);
599          } else {
600            _flow[e] -= d;
601            _res_cap[e] += d;
602            u = _graph.target(e);
603          }
604        }
605        _excess[s] -= d;
606        _excess[t] += d;
607      }
608
609      // Handling non-zero lower bounds
610      if (_lower) {
611        for (EdgeIt e(_graph); e != INVALID; ++e)
612          _flow[e] += (*_lower)[e];
613      }
614      return true;
615    }
616
617  }; //class CapacityScaling
618
619  ///@}
620
621} //namespace lemon
622
623#endif //LEMON_CAPACITY_SCALING_H
Note: See TracBrowser for help on using the repository browser.