COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/capacity_scaling.h @ 2581:054566ac0934

Last change on this file since 2581:054566ac0934 was 2581:054566ac0934, checked in by Peter Kovacs, 11 years ago

Query improvements in the min cost flow algorithms.

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