gravatar
alpar (Alpar Juttner)
alpar@cs.elte.hu
Merge #181, #323
0 2 0
merge default
0 files changed with 217 insertions and 95 deletions:
↑ Collapse diff ↑
Ignore white space 12 line context
... ...
@@ -26,12 +26,13 @@
26 26

	
27 27
#include <vector>
28 28
#include <limits>
29 29
#include <lemon/bin_heap.h>
30 30
#include <lemon/path.h>
31 31
#include <lemon/list_graph.h>
32
#include <lemon/dijkstra.h>
32 33
#include <lemon/maps.h>
33 34

	
34 35
namespace lemon {
35 36

	
36 37
  /// \addtogroup shortest_path
37 38
  /// @{
... ...
@@ -43,24 +44,24 @@
43 44
  /// finding arc-disjoint paths having minimum total length (cost)
44 45
  /// from a given source node to a given target node in a digraph.
45 46
  ///
46 47
  /// Note that this problem is a special case of the \ref min_cost_flow
47 48
  /// "minimum cost flow problem". This implementation is actually an
48 49
  /// efficient specialized version of the \ref CapacityScaling
49
  /// "Successive Shortest Path" algorithm directly for this problem.
50
  /// "successive shortest path" algorithm directly for this problem.
50 51
  /// Therefore this class provides query functions for flow values and
51 52
  /// node potentials (the dual solution) just like the minimum cost flow
52 53
  /// algorithms.
53 54
  ///
54 55
  /// \tparam GR The digraph type the algorithm runs on.
55 56
  /// \tparam LEN The type of the length map.
56 57
  /// The default value is <tt>GR::ArcMap<int></tt>.
57 58
  ///
58
  /// \warning Length values should be \e non-negative \e integers.
59
  /// \warning Length values should be \e non-negative.
59 60
  ///
60
  /// \note For finding node-disjoint paths this algorithm can be used
61
  /// \note For finding \e node-disjoint paths, this algorithm can be used
61 62
  /// along with the \ref SplitNodes adaptor.
62 63
#ifdef DOXYGEN
63 64
  template <typename GR, typename LEN>
64 65
#else
65 66
  template < typename GR,
66 67
             typename LEN = typename GR::template ArcMap<int> >
... ...
@@ -94,121 +95,163 @@
94 95

	
95 96
    /// The type of the path structures.
96 97
    typedef SimplePath<GR> Path;
97 98

	
98 99
  private:
99 100

	
101
    typedef typename Digraph::template NodeMap<int> HeapCrossRef;
102
    typedef BinHeap<Length, HeapCrossRef> Heap;
103

	
100 104
    // ResidualDijkstra is a special implementation of the
101 105
    // Dijkstra algorithm for finding shortest paths in the
102 106
    // residual network with respect to the reduced arc lengths
103 107
    // and modifying the node potentials according to the
104 108
    // distance of the nodes.
105 109
    class ResidualDijkstra
106 110
    {
107
      typedef typename Digraph::template NodeMap<int> HeapCrossRef;
108
      typedef BinHeap<Length, HeapCrossRef> Heap;
109

	
110 111
    private:
111 112

	
112
      // The digraph the algorithm runs on
113 113
      const Digraph &_graph;
114

	
115
      // The main maps
114
      const LengthMap &_length;
116 115
      const FlowMap &_flow;
117
      const LengthMap &_length;
118
      PotentialMap &_potential;
119

	
120
      // The distance map
121
      PotentialMap _dist;
122
      // The pred arc map
116
      PotentialMap &_pi;
123 117
      PredMap &_pred;
124
      // The processed (i.e. permanently labeled) nodes
125
      std::vector<Node> _proc_nodes;
126

	
127 118
      Node _s;
128 119
      Node _t;
120
      
121
      PotentialMap _dist;
122
      std::vector<Node> _proc_nodes;
129 123

	
130 124
    public:
131 125

	
132
      /// Constructor.
133
      ResidualDijkstra( const Digraph &graph,
134
                        const FlowMap &flow,
135
                        const LengthMap &length,
136
                        PotentialMap &potential,
137
                        PredMap &pred,
138
                        Node s, Node t ) :
139
        _graph(graph), _flow(flow), _length(length), _potential(potential),
140
        _dist(graph), _pred(pred), _s(s), _t(t) {}
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
      }
141 137

	
142
      /// \brief Run the algorithm. It returns \c true if a path is found
143
      /// from the source node to the target node.
144
      bool run() {
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() {
145 143
        HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
146 144
        Heap heap(heap_cross_ref);
147 145
        heap.push(_s, 0);
148 146
        _pred[_s] = INVALID;
149 147
        _proc_nodes.clear();
150 148

	
151 149
        // Process nodes
152 150
        while (!heap.empty() && heap.top() != _t) {
153 151
          Node u = heap.top(), v;
154
          Length d = heap.prio() + _potential[u], nd;
152
          Length d = heap.prio(), dn;
155 153
          _dist[u] = heap.prio();
154
          _proc_nodes.push_back(u);
156 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();
157 199
          _proc_nodes.push_back(u);
200
          heap.pop();
158 201

	
159 202
          // Traverse outgoing arcs
160 203
          for (OutArcIt e(_graph, u); e != INVALID; ++e) {
161 204
            if (_flow[e] == 0) {
162 205
              v = _graph.target(e);
163 206
              switch(heap.state(v)) {
164
              case Heap::PRE_HEAP:
165
                heap.push(v, d + _length[e] - _potential[v]);
166
                _pred[v] = e;
167
                break;
168
              case Heap::IN_HEAP:
169
                nd = d + _length[e] - _potential[v];
170
                if (nd < heap[v]) {
171
                  heap.decrease(v, nd);
207
                case Heap::PRE_HEAP:
208
                  heap.push(v, d + _length[e] - _pi[v]);
172 209
                  _pred[v] = e;
173
                }
174
                break;
175
              case Heap::POST_HEAP:
176
                break;
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;
177 220
              }
178 221
            }
179 222
          }
180 223

	
181 224
          // Traverse incoming arcs
182 225
          for (InArcIt e(_graph, u); e != INVALID; ++e) {
183 226
            if (_flow[e] == 1) {
184 227
              v = _graph.source(e);
185 228
              switch(heap.state(v)) {
186
              case Heap::PRE_HEAP:
187
                heap.push(v, d - _length[e] - _potential[v]);
188
                _pred[v] = e;
189
                break;
190
              case Heap::IN_HEAP:
191
                nd = d - _length[e] - _potential[v];
192
                if (nd < heap[v]) {
193
                  heap.decrease(v, nd);
229
                case Heap::PRE_HEAP:
230
                  heap.push(v, d - _length[e] - _pi[v]);
194 231
                  _pred[v] = e;
195
                }
196
                break;
197
              case Heap::POST_HEAP:
198
                break;
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;
199 242
              }
200 243
            }
201 244
          }
202 245
        }
203 246
        if (heap.empty()) return false;
204 247

	
205 248
        // Update potentials of processed nodes
206 249
        Length t_dist = heap.prio();
207 250
        for (int i = 0; i < int(_proc_nodes.size()); ++i)
208
          _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
251
          _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
209 252
        return true;
210 253
      }
211 254

	
212 255
    }; //class ResidualDijkstra
213 256

	
214 257
  private:
... ...
@@ -223,48 +266,49 @@
223 266
    bool _local_flow;
224 267
    // Node map of the current potentials
225 268
    PotentialMap *_potential;
226 269
    bool _local_potential;
227 270

	
228 271
    // The source node
229
    Node _source;
272
    Node _s;
230 273
    // The target node
231
    Node _target;
274
    Node _t;
232 275

	
233 276
    // Container to store the found paths
234
    std::vector< SimplePath<Digraph> > paths;
277
    std::vector<Path> _paths;
235 278
    int _path_num;
236 279

	
237 280
    // The pred arc map
238 281
    PredMap _pred;
239
    // Implementation of the Dijkstra algorithm for finding augmenting
240
    // shortest paths in the residual network
241
    ResidualDijkstra *_dijkstra;
282
    
283
    // Data for full init
284
    PotentialMap *_init_dist;
285
    PredMap *_init_pred;
286
    bool _full_init;
242 287

	
243 288
  public:
244 289

	
245 290
    /// \brief Constructor.
246 291
    ///
247 292
    /// Constructor.
248 293
    ///
249 294
    /// \param graph The digraph the algorithm runs on.
250 295
    /// \param length The length (cost) values of the arcs.
251 296
    Suurballe( const Digraph &graph,
252 297
               const LengthMap &length ) :
253 298
      _graph(graph), _length(length), _flow(0), _local_flow(false),
254
      _potential(0), _local_potential(false), _pred(graph)
255
    {
256
      LEMON_ASSERT(std::numeric_limits<Length>::is_integer,
257
        "The length type of Suurballe must be integer");
258
    }
299
      _potential(0), _local_potential(false), _pred(graph),
300
      _init_dist(0), _init_pred(0)
301
    {}
259 302

	
260 303
    /// Destructor.
261 304
    ~Suurballe() {
262 305
      if (_local_flow) delete _flow;
263 306
      if (_local_potential) delete _potential;
264
      delete _dijkstra;
307
      delete _init_dist;
308
      delete _init_pred;
265 309
    }
266 310

	
267 311
    /// \brief Set the flow map.
268 312
    ///
269 313
    /// This function sets the flow map.
270 314
    /// If it is not used before calling \ref run() or \ref init(),
... ...
@@ -303,16 +347,19 @@
303 347
      _potential = &map;
304 348
      return *this;
305 349
    }
306 350

	
307 351
    /// \name Execution Control
308 352
    /// The simplest way to execute the algorithm is to call the run()
309
    /// function.
310
    /// \n
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
311 357
    /// If you only need the flow that is the union of the found
312
    /// arc-disjoint paths, you may call init() and findFlow().
358
    /// arc-disjoint paths, then you may call findFlow() instead of
359
    /// start().
313 360

	
314 361
    /// @{
315 362

	
316 363
    /// \brief Run the algorithm.
317 364
    ///
318 365
    /// This function runs the algorithm.
... ...
@@ -326,41 +373,94 @@
326 373
    /// arc-disjoint paths found.
327 374
    ///
328 375
    /// \note Apart from the return value, <tt>s.run(s, t, k)</tt> is
329 376
    /// just a shortcut of the following code.
330 377
    /// \code
331 378
    ///   s.init(s);
332
    ///   s.findFlow(t, k);
333
    ///   s.findPaths();
379
    ///   s.start(t, k);
334 380
    /// \endcode
335 381
    int run(const Node& s, const Node& t, int k = 2) {
336 382
      init(s);
337
      findFlow(t, k);
338
      findPaths();
383
      start(t, k);
339 384
      return _path_num;
340 385
    }
341 386

	
342 387
    /// \brief Initialize the algorithm.
343 388
    ///
344
    /// This function initializes the algorithm.
389
    /// This function initializes the algorithm with the given source node.
345 390
    ///
346 391
    /// \param s The source node.
347 392
    void init(const Node& s) {
348
      _source = s;
393
      _s = s;
349 394

	
350 395
      // Initialize maps
351 396
      if (!_flow) {
352 397
        _flow = new FlowMap(_graph);
353 398
        _local_flow = true;
354 399
      }
355 400
      if (!_potential) {
356 401
        _potential = new PotentialMap(_graph);
357 402
        _local_potential = true;
358 403
      }
359
      for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
360
      for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
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;
361 461
    }
362 462

	
363 463
    /// \brief Execute the algorithm to find an optimal flow.
364 464
    ///
365 465
    /// This function executes the successive shortest path algorithm to
366 466
    /// find a minimum cost flow, which is the union of \c k (or less)
... ...
@@ -372,26 +472,45 @@
372 472
    /// \return \c k if there are at least \c k arc-disjoint paths from
373 473
    /// the source node to the given node \c t in the digraph.
374 474
    /// Otherwise it returns the number of arc-disjoint paths found.
375 475
    ///
376 476
    /// \pre \ref init() must be called before using this function.
377 477
    int findFlow(const Node& t, int k = 2) {
378
      _target = t;
379
      _dijkstra =
380
        new ResidualDijkstra( _graph, *_flow, _length, *_potential, _pred,
381
                              _source, _target );
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
      }
382 502

	
383 503
      // Find shortest paths
384
      _path_num = 0;
385 504
      while (_path_num < k) {
386 505
        // Run Dijkstra
387
        if (!_dijkstra->run()) break;
506
        if (!dijkstra.run(_path_num)) break;
388 507
        ++_path_num;
389 508

	
390 509
        // Set the flow along the found shortest path
391
        Node u = _target;
510
        Node u = _t;
392 511
        Arc e;
393 512
        while ((e = _pred[u]) != INVALID) {
394 513
          if (u == _graph.target(e)) {
395 514
            (*_flow)[e] = 1;
396 515
            u = _graph.source(e);
397 516
          } else {
... ...
@@ -402,30 +521,30 @@
402 521
      }
403 522
      return _path_num;
404 523
    }
405 524

	
406 525
    /// \brief Compute the paths from the flow.
407 526
    ///
408
    /// This function computes the paths from the found minimum cost flow,
409
    /// which is the union of some arc-disjoint paths.
527
    /// This function computes arc-disjoint paths from the found minimum
528
    /// cost flow, which is the union of them.
410 529
    ///
411 530
    /// \pre \ref init() and \ref findFlow() must be called before using
412 531
    /// this function.
413 532
    void findPaths() {
414 533
      FlowMap res_flow(_graph);
415 534
      for(ArcIt a(_graph); a != INVALID; ++a) res_flow[a] = (*_flow)[a];
416 535

	
417
      paths.clear();
418
      paths.resize(_path_num);
536
      _paths.clear();
537
      _paths.resize(_path_num);
419 538
      for (int i = 0; i < _path_num; ++i) {
420
        Node n = _source;
421
        while (n != _target) {
539
        Node n = _s;
540
        while (n != _t) {
422 541
          OutArcIt e(_graph, n);
423 542
          for ( ; res_flow[e] == 0; ++e) ;
424 543
          n = _graph.target(e);
425
          paths[i].addBack(e);
544
          _paths[i].addBack(e);
426 545
          res_flow[e] = 0;
427 546
        }
428 547
      }
429 548
    }
430 549

	
431 550
    /// @}
... ...
@@ -517,14 +636,14 @@
517 636
    ///
518 637
    /// \param i The function returns the <tt>i</tt>-th path.
519 638
    /// \c i must be between \c 0 and <tt>%pathNum()-1</tt>.
520 639
    ///
521 640
    /// \pre \ref run() or \ref findPaths() must be called before using
522 641
    /// this function.
523
    Path path(int i) const {
524
      return paths[i];
642
    const Path& path(int i) const {
643
      return _paths[i];
525 644
    }
526 645

	
527 646
    /// @}
528 647

	
529 648
  }; //class Suurballe
530 649

	
Ignore white space 12 line context
... ...
@@ -98,12 +98,15 @@
98 98
    .potentialMap(pi);
99 99

	
100 100
  int k;
101 101
  k = suurb_test.run(n, n);
102 102
  k = suurb_test.run(n, n, k);
103 103
  suurb_test.init(n);
104
  suurb_test.fullInit(n);
105
  suurb_test.start(n);
106
  suurb_test.start(n, k);
104 107
  k = suurb_test.findFlow(n);
105 108
  k = suurb_test.findFlow(n, k);
106 109
  suurb_test.findPaths();
107 110
  
108 111
  int f;
109 112
  VType c;
0 comments (0 inline)