COIN-OR::LEMON - Graph Library

source: lemon/lemon/suurballe.h @ 927:9a7e4e606f83

Last change on this file since 927:9a7e4e606f83 was 927:9a7e4e606f83, checked in by Peter Kovacs <kpeter@…>, 14 years ago

Add a fullInit() function to Suurballe (#181, #323)
to provide faster handling of multiple targets.
A start() function is also added, just for the sake of
convenience.

File size: 19.4 KB
Line 
1/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library.
4 *
5 * Copyright (C) 2003-2009
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_SUURBALLE_H
20#define LEMON_SUURBALLE_H
21
22///\ingroup shortest_path
23///\file
24///\brief An algorithm for finding arc-disjoint paths between two
25/// nodes having minimum total length.
26
27#include <vector>
28#include <limits>
29#include <lemon/bin_heap.h>
30#include <lemon/path.h>
31#include <lemon/list_graph.h>
32#include <lemon/dijkstra.h>
33#include <lemon/maps.h>
34
35namespace lemon {
36
37  /// \addtogroup shortest_path
38  /// @{
39
40  /// \brief Algorithm for finding arc-disjoint paths between two nodes
41  /// having minimum total length.
42  ///
43  /// \ref lemon::Suurballe "Suurballe" implements an algorithm for
44  /// finding arc-disjoint paths having minimum total length (cost)
45  /// from a given source node to a given target node in a digraph.
46  ///
47  /// Note that this problem is a special case of the \ref min_cost_flow
48  /// "minimum cost flow problem". This implementation is actually an
49  /// efficient specialized version of the \ref CapacityScaling
50  /// "successive shortest path" algorithm directly for this problem.
51  /// Therefore this class provides query functions for flow values and
52  /// node potentials (the dual solution) just like the minimum cost flow
53  /// algorithms.
54  ///
55  /// \tparam GR The digraph type the algorithm runs on.
56  /// \tparam LEN The type of the length map.
57  /// The default value is <tt>GR::ArcMap<int></tt>.
58  ///
59  /// \warning Length values should be \e non-negative.
60  ///
61  /// \note For finding \e node-disjoint paths, this algorithm can be used
62  /// along with the \ref SplitNodes adaptor.
63#ifdef DOXYGEN
64  template <typename GR, typename LEN>
65#else
66  template < typename GR,
67             typename LEN = typename GR::template ArcMap<int> >
68#endif
69  class Suurballe
70  {
71    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
72
73    typedef ConstMap<Arc, int> ConstArcMap;
74    typedef typename GR::template NodeMap<Arc> PredMap;
75
76  public:
77
78    /// The type of the digraph the algorithm runs on.
79    typedef GR Digraph;
80    /// The type of the length map.
81    typedef LEN LengthMap;
82    /// The type of the lengths.
83    typedef typename LengthMap::Value Length;
84#ifdef DOXYGEN
85    /// The type of the flow map.
86    typedef GR::ArcMap<int> FlowMap;
87    /// The type of the potential map.
88    typedef GR::NodeMap<Length> PotentialMap;
89#else
90    /// The type of the flow map.
91    typedef typename Digraph::template ArcMap<int> FlowMap;
92    /// The type of the potential map.
93    typedef typename Digraph::template NodeMap<Length> PotentialMap;
94#endif
95
96    /// The type of the path structures.
97    typedef SimplePath<GR> Path;
98
99  private:
100
101    typedef typename Digraph::template NodeMap<int> HeapCrossRef;
102    typedef BinHeap<Length, HeapCrossRef> Heap;
103
104    // ResidualDijkstra is a special implementation of the
105    // Dijkstra algorithm for finding shortest paths in the
106    // residual network with respect to the reduced arc lengths
107    // and modifying the node potentials according to the
108    // distance of the nodes.
109    class ResidualDijkstra
110    {
111    private:
112
113      const Digraph &_graph;
114      const LengthMap &_length;
115      const FlowMap &_flow;
116      PotentialMap &_pi;
117      PredMap &_pred;
118      Node _s;
119      Node _t;
120     
121      PotentialMap _dist;
122      std::vector<Node> _proc_nodes;
123
124    public:
125
126      // Constructor
127      ResidualDijkstra(Suurballe &srb) :
128        _graph(srb._graph), _length(srb._length),
129        _flow(*srb._flow), _pi(*srb._potential), _pred(srb._pred),
130        _s(srb._s), _t(srb._t), _dist(_graph) {}
131       
132      // Run the algorithm and return true if a path is found
133      // from the source node to the target node.
134      bool run(int cnt) {
135        return cnt == 0 ? startFirst() : start();
136      }
137
138    private:
139   
140      // Execute the algorithm for the first time (the flow and potential
141      // functions have to be identically zero).
142      bool startFirst() {
143        HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
144        Heap heap(heap_cross_ref);
145        heap.push(_s, 0);
146        _pred[_s] = INVALID;
147        _proc_nodes.clear();
148
149        // Process nodes
150        while (!heap.empty() && heap.top() != _t) {
151          Node u = heap.top(), v;
152          Length d = heap.prio(), dn;
153          _dist[u] = heap.prio();
154          _proc_nodes.push_back(u);
155          heap.pop();
156
157          // Traverse outgoing arcs
158          for (OutArcIt e(_graph, u); e != INVALID; ++e) {
159            v = _graph.target(e);
160            switch(heap.state(v)) {
161              case Heap::PRE_HEAP:
162                heap.push(v, d + _length[e]);
163                _pred[v] = e;
164                break;
165              case Heap::IN_HEAP:
166                dn = d + _length[e];
167                if (dn < heap[v]) {
168                  heap.decrease(v, dn);
169                  _pred[v] = e;
170                }
171                break;
172              case Heap::POST_HEAP:
173                break;
174            }
175          }
176        }
177        if (heap.empty()) return false;
178
179        // Update potentials of processed nodes
180        Length t_dist = heap.prio();
181        for (int i = 0; i < int(_proc_nodes.size()); ++i)
182          _pi[_proc_nodes[i]] = _dist[_proc_nodes[i]] - t_dist;
183        return true;
184      }
185
186      // Execute the algorithm.
187      bool start() {
188        HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
189        Heap heap(heap_cross_ref);
190        heap.push(_s, 0);
191        _pred[_s] = INVALID;
192        _proc_nodes.clear();
193
194        // Process nodes
195        while (!heap.empty() && heap.top() != _t) {
196          Node u = heap.top(), v;
197          Length d = heap.prio() + _pi[u], dn;
198          _dist[u] = heap.prio();
199          _proc_nodes.push_back(u);
200          heap.pop();
201
202          // Traverse outgoing arcs
203          for (OutArcIt e(_graph, u); e != INVALID; ++e) {
204            if (_flow[e] == 0) {
205              v = _graph.target(e);
206              switch(heap.state(v)) {
207                case Heap::PRE_HEAP:
208                  heap.push(v, d + _length[e] - _pi[v]);
209                  _pred[v] = e;
210                  break;
211                case Heap::IN_HEAP:
212                  dn = d + _length[e] - _pi[v];
213                  if (dn < heap[v]) {
214                    heap.decrease(v, dn);
215                    _pred[v] = e;
216                  }
217                  break;
218                case Heap::POST_HEAP:
219                  break;
220              }
221            }
222          }
223
224          // Traverse incoming arcs
225          for (InArcIt e(_graph, u); e != INVALID; ++e) {
226            if (_flow[e] == 1) {
227              v = _graph.source(e);
228              switch(heap.state(v)) {
229                case Heap::PRE_HEAP:
230                  heap.push(v, d - _length[e] - _pi[v]);
231                  _pred[v] = e;
232                  break;
233                case Heap::IN_HEAP:
234                  dn = d - _length[e] - _pi[v];
235                  if (dn < heap[v]) {
236                    heap.decrease(v, dn);
237                    _pred[v] = e;
238                  }
239                  break;
240                case Heap::POST_HEAP:
241                  break;
242              }
243            }
244          }
245        }
246        if (heap.empty()) return false;
247
248        // Update potentials of processed nodes
249        Length t_dist = heap.prio();
250        for (int i = 0; i < int(_proc_nodes.size()); ++i)
251          _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
252        return true;
253      }
254
255    }; //class ResidualDijkstra
256
257  private:
258
259    // The digraph the algorithm runs on
260    const Digraph &_graph;
261    // The length map
262    const LengthMap &_length;
263
264    // Arc map of the current flow
265    FlowMap *_flow;
266    bool _local_flow;
267    // Node map of the current potentials
268    PotentialMap *_potential;
269    bool _local_potential;
270
271    // The source node
272    Node _s;
273    // The target node
274    Node _t;
275
276    // Container to store the found paths
277    std::vector<Path> _paths;
278    int _path_num;
279
280    // The pred arc map
281    PredMap _pred;
282   
283    // Data for full init
284    PotentialMap *_init_dist;
285    PredMap *_init_pred;
286    bool _full_init;
287
288  public:
289
290    /// \brief Constructor.
291    ///
292    /// Constructor.
293    ///
294    /// \param graph The digraph the algorithm runs on.
295    /// \param length The length (cost) values of the arcs.
296    Suurballe( const Digraph &graph,
297               const LengthMap &length ) :
298      _graph(graph), _length(length), _flow(0), _local_flow(false),
299      _potential(0), _local_potential(false), _pred(graph),
300      _init_dist(0), _init_pred(0)
301    {}
302
303    /// Destructor.
304    ~Suurballe() {
305      if (_local_flow) delete _flow;
306      if (_local_potential) delete _potential;
307      delete _init_dist;
308      delete _init_pred;
309    }
310
311    /// \brief Set the flow map.
312    ///
313    /// This function sets the flow map.
314    /// If it is not used before calling \ref run() or \ref init(),
315    /// an instance will be allocated automatically. The destructor
316    /// deallocates this automatically allocated map, of course.
317    ///
318    /// The found flow contains only 0 and 1 values, since it is the
319    /// union of the found arc-disjoint paths.
320    ///
321    /// \return <tt>(*this)</tt>
322    Suurballe& flowMap(FlowMap &map) {
323      if (_local_flow) {
324        delete _flow;
325        _local_flow = false;
326      }
327      _flow = &map;
328      return *this;
329    }
330
331    /// \brief Set the potential map.
332    ///
333    /// This function sets the potential map.
334    /// If it is not used before calling \ref run() or \ref init(),
335    /// an instance will be allocated automatically. The destructor
336    /// deallocates this automatically allocated map, of course.
337    ///
338    /// The node potentials provide the dual solution of the underlying
339    /// \ref min_cost_flow "minimum cost flow problem".
340    ///
341    /// \return <tt>(*this)</tt>
342    Suurballe& potentialMap(PotentialMap &map) {
343      if (_local_potential) {
344        delete _potential;
345        _local_potential = false;
346      }
347      _potential = &map;
348      return *this;
349    }
350
351    /// \name Execution Control
352    /// The simplest way to execute the algorithm is to call the run()
353    /// function.\n
354    /// If you need to execute the algorithm many times using the same
355    /// source node, then you may call fullInit() once and start()
356    /// for each target node.\n
357    /// If you only need the flow that is the union of the found
358    /// arc-disjoint paths, then you may call findFlow() instead of
359    /// start().
360
361    /// @{
362
363    /// \brief Run the algorithm.
364    ///
365    /// This function runs the algorithm.
366    ///
367    /// \param s The source node.
368    /// \param t The target node.
369    /// \param k The number of paths to be found.
370    ///
371    /// \return \c k if there are at least \c k arc-disjoint paths from
372    /// \c s to \c t in the digraph. Otherwise it returns the number of
373    /// arc-disjoint paths found.
374    ///
375    /// \note Apart from the return value, <tt>s.run(s, t, k)</tt> is
376    /// just a shortcut of the following code.
377    /// \code
378    ///   s.init(s);
379    ///   s.start(t, k);
380    /// \endcode
381    int run(const Node& s, const Node& t, int k = 2) {
382      init(s);
383      start(t, k);
384      return _path_num;
385    }
386
387    /// \brief Initialize the algorithm.
388    ///
389    /// This function initializes the algorithm with the given source node.
390    ///
391    /// \param s The source node.
392    void init(const Node& s) {
393      _s = s;
394
395      // Initialize maps
396      if (!_flow) {
397        _flow = new FlowMap(_graph);
398        _local_flow = true;
399      }
400      if (!_potential) {
401        _potential = new PotentialMap(_graph);
402        _local_potential = true;
403      }
404      _full_init = false;
405    }
406
407    /// \brief Initialize the algorithm and perform Dijkstra.
408    ///
409    /// This function initializes the algorithm and performs a full
410    /// Dijkstra search from the given source node. It makes consecutive
411    /// executions of \ref start() "start(t, k)" faster, since they
412    /// have to perform %Dijkstra only k-1 times.
413    ///
414    /// This initialization is usually worth using instead of \ref init()
415    /// if the algorithm is executed many times using the same source node.
416    ///
417    /// \param s The source node.
418    void fullInit(const Node& s) {
419      // Initialize maps
420      init(s);
421      if (!_init_dist) {
422        _init_dist = new PotentialMap(_graph);
423      }
424      if (!_init_pred) {
425        _init_pred = new PredMap(_graph);
426      }
427
428      // Run a full Dijkstra
429      typename Dijkstra<Digraph, LengthMap>
430        ::template SetStandardHeap<Heap>
431        ::template SetDistMap<PotentialMap>
432        ::template SetPredMap<PredMap>
433        ::Create dijk(_graph, _length);
434      dijk.distMap(*_init_dist).predMap(*_init_pred);
435      dijk.run(s);
436     
437      _full_init = true;
438    }
439
440    /// \brief Execute the algorithm.
441    ///
442    /// This function executes the algorithm.
443    ///
444    /// \param t The target node.
445    /// \param k The number of paths to be found.
446    ///
447    /// \return \c k if there are at least \c k arc-disjoint paths from
448    /// \c s to \c t in the digraph. Otherwise it returns the number of
449    /// arc-disjoint paths found.
450    ///
451    /// \note Apart from the return value, <tt>s.start(t, k)</tt> is
452    /// just a shortcut of the following code.
453    /// \code
454    ///   s.findFlow(t, k);
455    ///   s.findPaths();
456    /// \endcode
457    int start(const Node& t, int k = 2) {
458      findFlow(t, k);
459      findPaths();
460      return _path_num;
461    }
462
463    /// \brief Execute the algorithm to find an optimal flow.
464    ///
465    /// This function executes the successive shortest path algorithm to
466    /// find a minimum cost flow, which is the union of \c k (or less)
467    /// arc-disjoint paths.
468    ///
469    /// \param t The target node.
470    /// \param k The number of paths to be found.
471    ///
472    /// \return \c k if there are at least \c k arc-disjoint paths from
473    /// the source node to the given node \c t in the digraph.
474    /// Otherwise it returns the number of arc-disjoint paths found.
475    ///
476    /// \pre \ref init() must be called before using this function.
477    int findFlow(const Node& t, int k = 2) {
478      _t = t;
479      ResidualDijkstra dijkstra(*this);
480     
481      // Initialization
482      for (ArcIt e(_graph); e != INVALID; ++e) {
483        (*_flow)[e] = 0;
484      }
485      if (_full_init) {
486        for (NodeIt n(_graph); n != INVALID; ++n) {
487          (*_potential)[n] = (*_init_dist)[n];
488        }
489        Node u = _t;
490        Arc e;
491        while ((e = (*_init_pred)[u]) != INVALID) {
492          (*_flow)[e] = 1;
493          u = _graph.source(e);
494        }       
495        _path_num = 1;
496      } else {
497        for (NodeIt n(_graph); n != INVALID; ++n) {
498          (*_potential)[n] = 0;
499        }
500        _path_num = 0;
501      }
502
503      // Find shortest paths
504      while (_path_num < k) {
505        // Run Dijkstra
506        if (!dijkstra.run(_path_num)) break;
507        ++_path_num;
508
509        // Set the flow along the found shortest path
510        Node u = _t;
511        Arc e;
512        while ((e = _pred[u]) != INVALID) {
513          if (u == _graph.target(e)) {
514            (*_flow)[e] = 1;
515            u = _graph.source(e);
516          } else {
517            (*_flow)[e] = 0;
518            u = _graph.target(e);
519          }
520        }
521      }
522      return _path_num;
523    }
524
525    /// \brief Compute the paths from the flow.
526    ///
527    /// This function computes arc-disjoint paths from the found minimum
528    /// cost flow, which is the union of them.
529    ///
530    /// \pre \ref init() and \ref findFlow() must be called before using
531    /// this function.
532    void findPaths() {
533      FlowMap res_flow(_graph);
534      for(ArcIt a(_graph); a != INVALID; ++a) res_flow[a] = (*_flow)[a];
535
536      _paths.clear();
537      _paths.resize(_path_num);
538      for (int i = 0; i < _path_num; ++i) {
539        Node n = _s;
540        while (n != _t) {
541          OutArcIt e(_graph, n);
542          for ( ; res_flow[e] == 0; ++e) ;
543          n = _graph.target(e);
544          _paths[i].addBack(e);
545          res_flow[e] = 0;
546        }
547      }
548    }
549
550    /// @}
551
552    /// \name Query Functions
553    /// The results of the algorithm can be obtained using these
554    /// functions.
555    /// \n The algorithm should be executed before using them.
556
557    /// @{
558
559    /// \brief Return the total length of the found paths.
560    ///
561    /// This function returns the total length of the found paths, i.e.
562    /// the total cost of the found flow.
563    /// The complexity of the function is O(e).
564    ///
565    /// \pre \ref run() or \ref findFlow() must be called before using
566    /// this function.
567    Length totalLength() const {
568      Length c = 0;
569      for (ArcIt e(_graph); e != INVALID; ++e)
570        c += (*_flow)[e] * _length[e];
571      return c;
572    }
573
574    /// \brief Return the flow value on the given arc.
575    ///
576    /// This function returns the flow value on the given arc.
577    /// It is \c 1 if the arc is involved in one of the found arc-disjoint
578    /// paths, otherwise it is \c 0.
579    ///
580    /// \pre \ref run() or \ref findFlow() must be called before using
581    /// this function.
582    int flow(const Arc& arc) const {
583      return (*_flow)[arc];
584    }
585
586    /// \brief Return a const reference to an arc map storing the
587    /// found flow.
588    ///
589    /// This function returns a const reference to an arc map storing
590    /// the flow that is the union of the found arc-disjoint paths.
591    ///
592    /// \pre \ref run() or \ref findFlow() must be called before using
593    /// this function.
594    const FlowMap& flowMap() const {
595      return *_flow;
596    }
597
598    /// \brief Return the potential of the given node.
599    ///
600    /// This function returns the potential of the given node.
601    /// The node potentials provide the dual solution of the
602    /// underlying \ref min_cost_flow "minimum cost flow problem".
603    ///
604    /// \pre \ref run() or \ref findFlow() must be called before using
605    /// this function.
606    Length potential(const Node& node) const {
607      return (*_potential)[node];
608    }
609
610    /// \brief Return a const reference to a node map storing the
611    /// found potentials (the dual solution).
612    ///
613    /// This function returns a const reference to a node map storing
614    /// the found potentials that provide the dual solution of the
615    /// underlying \ref min_cost_flow "minimum cost flow problem".
616    ///
617    /// \pre \ref run() or \ref findFlow() must be called before using
618    /// this function.
619    const PotentialMap& potentialMap() const {
620      return *_potential;
621    }
622
623    /// \brief Return the number of the found paths.
624    ///
625    /// This function returns the number of the found paths.
626    ///
627    /// \pre \ref run() or \ref findFlow() must be called before using
628    /// this function.
629    int pathNum() const {
630      return _path_num;
631    }
632
633    /// \brief Return a const reference to the specified path.
634    ///
635    /// This function returns a const reference to the specified path.
636    ///
637    /// \param i The function returns the <tt>i</tt>-th path.
638    /// \c i must be between \c 0 and <tt>%pathNum()-1</tt>.
639    ///
640    /// \pre \ref run() or \ref findPaths() must be called before using
641    /// this function.
642    const Path& path(int i) const {
643      return _paths[i];
644    }
645
646    /// @}
647
648  }; //class Suurballe
649
650  ///@}
651
652} //namespace lemon
653
654#endif //LEMON_SUURBALLE_H
Note: See TracBrowser for help on using the repository browser.