COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/capacity_scaling.h @ 2535:716024e7c080

Last change on this file since 2535:716024e7c080 was 2535:716024e7c080, checked in by Peter Kovacs, 11 years ago

Redesigned CapacityScaling? algorithm with almost the same interface.
The new version does not use the ResidualGraphAdaptor? for performance reasons.
Scaling can be enabled and disabled with a parameter of the run() function.

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