COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/capacity_scaling.h @ 2589:1bbb28acb8c9

Last change on this file since 2589:1bbb28acb8c9 was 2589:1bbb28acb8c9, checked in by Peter Kovacs, 12 years ago

Bug fix in CapacityScaling?.

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