COIN-OR::LEMON - Graph Library

source: lemon/lemon/capacity_scaling.h @ 871:d3e32a777d0b

Last change on this file since 871:d3e32a777d0b was 871:d3e32a777d0b, checked in by Peter Kovacs <kpeter@…>, 10 years ago

Port CapacityScaling? from SVN -r3524 (#180)

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