gravatar
deba@inf.elte.hu
deba@inf.elte.hu
Fix typos (ticket #192)
0 3 0
default
3 files changed with 4 insertions and 4 deletions:
↑ Collapse diff ↑
Ignore white space 1536 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2008
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_CONNECTIVITY_H
20 20
#define LEMON_CONNECTIVITY_H
21 21

	
22 22
#include <lemon/dfs.h>
23 23
#include <lemon/bfs.h>
24 24
#include <lemon/core.h>
25 25
#include <lemon/maps.h>
26 26
#include <lemon/adaptors.h>
27 27

	
28 28
#include <lemon/concepts/digraph.h>
29 29
#include <lemon/concepts/graph.h>
30 30
#include <lemon/concept_check.h>
31 31

	
32 32
#include <stack>
33 33
#include <functional>
34 34

	
35 35
/// \ingroup connectivity
36 36
/// \file
37 37
/// \brief Connectivity algorithms
38 38
///
39 39
/// Connectivity algorithms
40 40

	
41 41
namespace lemon {
42 42

	
43 43
  /// \ingroup connectivity
44 44
  ///
45 45
  /// \brief Check whether the given undirected graph is connected.
46 46
  ///
47 47
  /// Check whether the given undirected graph is connected.
48 48
  /// \param graph The undirected graph.
49 49
  /// \return %True when there is path between any two nodes in the graph.
50 50
  /// \note By definition, the empty graph is connected.
51 51
  template <typename Graph>
52 52
  bool connected(const Graph& graph) {
53 53
    checkConcept<concepts::Graph, Graph>();
54 54
    typedef typename Graph::NodeIt NodeIt;
55 55
    if (NodeIt(graph) == INVALID) return true;
56 56
    Dfs<Graph> dfs(graph);
57 57
    dfs.run(NodeIt(graph));
58 58
    for (NodeIt it(graph); it != INVALID; ++it) {
59 59
      if (!dfs.reached(it)) {
60 60
        return false;
61 61
      }
62 62
    }
63 63
    return true;
64 64
  }
65 65

	
66 66
  /// \ingroup connectivity
67 67
  ///
68 68
  /// \brief Count the number of connected components of an undirected graph
69 69
  ///
70 70
  /// Count the number of connected components of an undirected graph
71 71
  ///
72 72
  /// \param graph The graph. It must be undirected.
73 73
  /// \return The number of components
74 74
  /// \note By definition, the empty graph consists
75 75
  /// of zero connected components.
76 76
  template <typename Graph>
77 77
  int countConnectedComponents(const Graph &graph) {
78 78
    checkConcept<concepts::Graph, Graph>();
79 79
    typedef typename Graph::Node Node;
80 80
    typedef typename Graph::Arc Arc;
81 81

	
82 82
    typedef NullMap<Node, Arc> PredMap;
83 83
    typedef NullMap<Node, int> DistMap;
84 84

	
85 85
    int compNum = 0;
86 86
    typename Bfs<Graph>::
87 87
      template SetPredMap<PredMap>::
88 88
      template SetDistMap<DistMap>::
89 89
      Create bfs(graph);
90 90

	
91 91
    PredMap predMap;
92 92
    bfs.predMap(predMap);
93 93

	
94 94
    DistMap distMap;
95 95
    bfs.distMap(distMap);
96 96

	
97 97
    bfs.init();
98 98
    for(typename Graph::NodeIt n(graph); n != INVALID; ++n) {
99 99
      if (!bfs.reached(n)) {
100 100
        bfs.addSource(n);
101 101
        bfs.start();
102 102
        ++compNum;
103 103
      }
104 104
    }
105 105
    return compNum;
106 106
  }
107 107

	
108 108
  /// \ingroup connectivity
109 109
  ///
110 110
  /// \brief Find the connected components of an undirected graph
111 111
  ///
112 112
  /// Find the connected components of an undirected graph.
113 113
  ///
114 114
  /// \param graph The graph. It must be undirected.
115 115
  /// \retval compMap A writable node map. The values will be set from 0 to
116 116
  /// the number of the connected components minus one. Each values of the map
117 117
  /// will be set exactly once, the values of a certain component will be
118 118
  /// set continuously.
119 119
  /// \return The number of components
120 120
  ///
121 121
  template <class Graph, class NodeMap>
122 122
  int connectedComponents(const Graph &graph, NodeMap &compMap) {
123 123
    checkConcept<concepts::Graph, Graph>();
124 124
    typedef typename Graph::Node Node;
125 125
    typedef typename Graph::Arc Arc;
126 126
    checkConcept<concepts::WriteMap<Node, int>, NodeMap>();
127 127

	
128 128
    typedef NullMap<Node, Arc> PredMap;
129 129
    typedef NullMap<Node, int> DistMap;
130 130

	
131 131
    int compNum = 0;
132 132
    typename Bfs<Graph>::
133 133
      template SetPredMap<PredMap>::
134 134
      template SetDistMap<DistMap>::
135 135
      Create bfs(graph);
136 136

	
137 137
    PredMap predMap;
138 138
    bfs.predMap(predMap);
139 139

	
140 140
    DistMap distMap;
141 141
    bfs.distMap(distMap);
142 142

	
143 143
    bfs.init();
144 144
    for(typename Graph::NodeIt n(graph); n != INVALID; ++n) {
145 145
      if(!bfs.reached(n)) {
146 146
        bfs.addSource(n);
147 147
        while (!bfs.emptyQueue()) {
148 148
          compMap.set(bfs.nextNode(), compNum);
149 149
          bfs.processNextNode();
150 150
        }
151 151
        ++compNum;
152 152
      }
153 153
    }
154 154
    return compNum;
155 155
  }
156 156

	
157 157
  namespace _connectivity_bits {
158 158

	
159 159
    template <typename Digraph, typename Iterator >
160 160
    struct LeaveOrderVisitor : public DfsVisitor<Digraph> {
161 161
    public:
162 162
      typedef typename Digraph::Node Node;
163 163
      LeaveOrderVisitor(Iterator it) : _it(it) {}
164 164

	
165 165
      void leave(const Node& node) {
166 166
        *(_it++) = node;
167 167
      }
168 168

	
169 169
    private:
170 170
      Iterator _it;
171 171
    };
172 172

	
173 173
    template <typename Digraph, typename Map>
174 174
    struct FillMapVisitor : public DfsVisitor<Digraph> {
175 175
    public:
176 176
      typedef typename Digraph::Node Node;
177 177
      typedef typename Map::Value Value;
178 178

	
179 179
      FillMapVisitor(Map& map, Value& value)
180 180
        : _map(map), _value(value) {}
181 181

	
182 182
      void reach(const Node& node) {
183 183
        _map.set(node, _value);
184 184
      }
185 185
    private:
186 186
      Map& _map;
187 187
      Value& _value;
188 188
    };
189 189

	
190 190
    template <typename Digraph, typename ArcMap>
191 191
    struct StronglyConnectedCutArcsVisitor : public DfsVisitor<Digraph> {
192 192
    public:
193 193
      typedef typename Digraph::Node Node;
194 194
      typedef typename Digraph::Arc Arc;
195 195

	
196 196
      StronglyConnectedCutArcsVisitor(const Digraph& digraph,
197 197
                                      ArcMap& cutMap,
198 198
                                      int& cutNum)
199 199
        : _digraph(digraph), _cutMap(cutMap), _cutNum(cutNum),
200 200
          _compMap(digraph, -1), _num(-1) {
201 201
      }
202 202

	
203 203
      void start(const Node&) {
204 204
        ++_num;
205 205
      }
206 206

	
207 207
      void reach(const Node& node) {
208 208
        _compMap.set(node, _num);
209 209
      }
210 210

	
211 211
      void examine(const Arc& arc) {
212 212
         if (_compMap[_digraph.source(arc)] !=
213 213
             _compMap[_digraph.target(arc)]) {
214 214
           _cutMap.set(arc, true);
215 215
           ++_cutNum;
216 216
         }
217 217
      }
218 218
    private:
219 219
      const Digraph& _digraph;
220 220
      ArcMap& _cutMap;
221 221
      int& _cutNum;
222 222

	
223 223
      typename Digraph::template NodeMap<int> _compMap;
224 224
      int _num;
225 225
    };
226 226

	
227 227
  }
228 228

	
229 229

	
230 230
  /// \ingroup connectivity
231 231
  ///
232 232
  /// \brief Check whether the given directed graph is strongly connected.
233 233
  ///
234 234
  /// Check whether the given directed graph is strongly connected. The
235 235
  /// graph is strongly connected when any two nodes of the graph are
236 236
  /// connected with directed paths in both direction.
237 237
  /// \return %False when the graph is not strongly connected.
238 238
  /// \see connected
239 239
  ///
240 240
  /// \note By definition, the empty graph is strongly connected.
241 241
  template <typename Digraph>
242 242
  bool stronglyConnected(const Digraph& digraph) {
243 243
    checkConcept<concepts::Digraph, Digraph>();
244 244

	
245 245
    typedef typename Digraph::Node Node;
246 246
    typedef typename Digraph::NodeIt NodeIt;
247 247

	
248 248
    typename Digraph::Node source = NodeIt(digraph);
249 249
    if (source == INVALID) return true;
250 250

	
251 251
    using namespace _connectivity_bits;
252 252

	
253 253
    typedef DfsVisitor<Digraph> Visitor;
254 254
    Visitor visitor;
255 255

	
256 256
    DfsVisit<Digraph, Visitor> dfs(digraph, visitor);
257 257
    dfs.init();
258 258
    dfs.addSource(source);
259 259
    dfs.start();
260 260

	
261 261
    for (NodeIt it(digraph); it != INVALID; ++it) {
262 262
      if (!dfs.reached(it)) {
263 263
        return false;
264 264
      }
265 265
    }
266 266

	
267 267
    typedef ReverseDigraph<const Digraph> RDigraph;
268 268
    typedef typename RDigraph::NodeIt RNodeIt;
269 269
    RDigraph rdigraph(digraph);
270 270

	
271 271
    typedef DfsVisitor<Digraph> RVisitor;
272 272
    RVisitor rvisitor;
273 273

	
274 274
    DfsVisit<RDigraph, RVisitor> rdfs(rdigraph, rvisitor);
275 275
    rdfs.init();
276 276
    rdfs.addSource(source);
277 277
    rdfs.start();
278 278

	
279 279
    for (RNodeIt it(rdigraph); it != INVALID; ++it) {
280 280
      if (!rdfs.reached(it)) {
281 281
        return false;
282 282
      }
283 283
    }
284 284

	
285 285
    return true;
286 286
  }
287 287

	
288 288
  /// \ingroup connectivity
289 289
  ///
290 290
  /// \brief Count the strongly connected components of a directed graph
291 291
  ///
292 292
  /// Count the strongly connected components of a directed graph.
293 293
  /// The strongly connected components are the classes of an
294 294
  /// equivalence relation on the nodes of the graph. Two nodes are in
295 295
  /// the same class if they are connected with directed paths in both
296 296
  /// direction.
297 297
  ///
298
  /// \param graph The graph.
298
  /// \param digraph The graph.
299 299
  /// \return The number of components
300 300
  /// \note By definition, the empty graph has zero
301 301
  /// strongly connected components.
302 302
  template <typename Digraph>
303 303
  int countStronglyConnectedComponents(const Digraph& digraph) {
304 304
    checkConcept<concepts::Digraph, Digraph>();
305 305

	
306 306
    using namespace _connectivity_bits;
307 307

	
308 308
    typedef typename Digraph::Node Node;
309 309
    typedef typename Digraph::Arc Arc;
310 310
    typedef typename Digraph::NodeIt NodeIt;
311 311
    typedef typename Digraph::ArcIt ArcIt;
312 312

	
313 313
    typedef std::vector<Node> Container;
314 314
    typedef typename Container::iterator Iterator;
315 315

	
316 316
    Container nodes(countNodes(digraph));
317 317
    typedef LeaveOrderVisitor<Digraph, Iterator> Visitor;
318 318
    Visitor visitor(nodes.begin());
319 319

	
320 320
    DfsVisit<Digraph, Visitor> dfs(digraph, visitor);
321 321
    dfs.init();
322 322
    for (NodeIt it(digraph); it != INVALID; ++it) {
323 323
      if (!dfs.reached(it)) {
324 324
        dfs.addSource(it);
325 325
        dfs.start();
326 326
      }
327 327
    }
328 328

	
329 329
    typedef typename Container::reverse_iterator RIterator;
330 330
    typedef ReverseDigraph<const Digraph> RDigraph;
331 331

	
332 332
    RDigraph rdigraph(digraph);
333 333

	
334 334
    typedef DfsVisitor<Digraph> RVisitor;
335 335
    RVisitor rvisitor;
336 336

	
337 337
    DfsVisit<RDigraph, RVisitor> rdfs(rdigraph, rvisitor);
338 338

	
339 339
    int compNum = 0;
340 340

	
341 341
    rdfs.init();
342 342
    for (RIterator it = nodes.rbegin(); it != nodes.rend(); ++it) {
343 343
      if (!rdfs.reached(*it)) {
344 344
        rdfs.addSource(*it);
345 345
        rdfs.start();
346 346
        ++compNum;
347 347
      }
348 348
    }
349 349
    return compNum;
350 350
  }
351 351

	
352 352
  /// \ingroup connectivity
353 353
  ///
354 354
  /// \brief Find the strongly connected components of a directed graph
355 355
  ///
356 356
  /// Find the strongly connected components of a directed graph.  The
357 357
  /// strongly connected components are the classes of an equivalence
358 358
  /// relation on the nodes of the graph. Two nodes are in
359 359
  /// relationship when there are directed paths between them in both
360 360
  /// direction. In addition, the numbering of components will satisfy
361 361
  /// that there is no arc going from a higher numbered component to
362 362
  /// a lower.
363 363
  ///
364 364
  /// \param digraph The digraph.
365 365
  /// \retval compMap A writable node map. The values will be set from 0 to
366 366
  /// the number of the strongly connected components minus one. Each value
367 367
  /// of the map will be set exactly once, the values of a certain component
368 368
  /// will be set continuously.
369 369
  /// \return The number of components
370 370
  ///
371 371
  template <typename Digraph, typename NodeMap>
372 372
  int stronglyConnectedComponents(const Digraph& digraph, NodeMap& compMap) {
373 373
    checkConcept<concepts::Digraph, Digraph>();
374 374
    typedef typename Digraph::Node Node;
375 375
    typedef typename Digraph::NodeIt NodeIt;
376 376
    checkConcept<concepts::WriteMap<Node, int>, NodeMap>();
377 377

	
378 378
    using namespace _connectivity_bits;
379 379

	
380 380
    typedef std::vector<Node> Container;
381 381
    typedef typename Container::iterator Iterator;
382 382

	
383 383
    Container nodes(countNodes(digraph));
384 384
    typedef LeaveOrderVisitor<Digraph, Iterator> Visitor;
385 385
    Visitor visitor(nodes.begin());
386 386

	
387 387
    DfsVisit<Digraph, Visitor> dfs(digraph, visitor);
388 388
    dfs.init();
389 389
    for (NodeIt it(digraph); it != INVALID; ++it) {
390 390
      if (!dfs.reached(it)) {
391 391
        dfs.addSource(it);
392 392
        dfs.start();
393 393
      }
394 394
    }
395 395

	
396 396
    typedef typename Container::reverse_iterator RIterator;
397 397
    typedef ReverseDigraph<const Digraph> RDigraph;
398 398

	
399 399
    RDigraph rdigraph(digraph);
400 400

	
401 401
    int compNum = 0;
402 402

	
403 403
    typedef FillMapVisitor<RDigraph, NodeMap> RVisitor;
404 404
    RVisitor rvisitor(compMap, compNum);
405 405

	
406 406
    DfsVisit<RDigraph, RVisitor> rdfs(rdigraph, rvisitor);
407 407

	
408 408
    rdfs.init();
409 409
    for (RIterator it = nodes.rbegin(); it != nodes.rend(); ++it) {
410 410
      if (!rdfs.reached(*it)) {
411 411
        rdfs.addSource(*it);
412 412
        rdfs.start();
413 413
        ++compNum;
414 414
      }
415 415
    }
416 416
    return compNum;
417 417
  }
418 418

	
419 419
  /// \ingroup connectivity
420 420
  ///
421 421
  /// \brief Find the cut arcs of the strongly connected components.
422 422
  ///
423 423
  /// Find the cut arcs of the strongly connected components.
424 424
  /// The strongly connected components are the classes of an equivalence
425 425
  /// relation on the nodes of the graph. Two nodes are in relationship
426 426
  /// when there are directed paths between them in both direction.
427 427
  /// The strongly connected components are separated by the cut arcs.
428 428
  ///
429 429
  /// \param graph The graph.
430 430
  /// \retval cutMap A writable node map. The values will be set true when the
431 431
  /// arc is a cut arc.
432 432
  ///
433 433
  /// \return The number of cut arcs
434 434
  template <typename Digraph, typename ArcMap>
435 435
  int stronglyConnectedCutArcs(const Digraph& graph, ArcMap& cutMap) {
436 436
    checkConcept<concepts::Digraph, Digraph>();
437 437
    typedef typename Digraph::Node Node;
438 438
    typedef typename Digraph::Arc Arc;
439 439
    typedef typename Digraph::NodeIt NodeIt;
440 440
    checkConcept<concepts::WriteMap<Arc, bool>, ArcMap>();
441 441

	
442 442
    using namespace _connectivity_bits;
443 443

	
444 444
    typedef std::vector<Node> Container;
445 445
    typedef typename Container::iterator Iterator;
446 446

	
447 447
    Container nodes(countNodes(graph));
448 448
    typedef LeaveOrderVisitor<Digraph, Iterator> Visitor;
449 449
    Visitor visitor(nodes.begin());
450 450

	
451 451
    DfsVisit<Digraph, Visitor> dfs(graph, visitor);
452 452
    dfs.init();
453 453
    for (NodeIt it(graph); it != INVALID; ++it) {
454 454
      if (!dfs.reached(it)) {
455 455
        dfs.addSource(it);
456 456
        dfs.start();
457 457
      }
458 458
    }
459 459

	
460 460
    typedef typename Container::reverse_iterator RIterator;
461 461
    typedef ReverseDigraph<const Digraph> RDigraph;
462 462

	
463 463
    RDigraph rgraph(graph);
464 464

	
465 465
    int cutNum = 0;
466 466

	
467 467
    typedef StronglyConnectedCutArcsVisitor<RDigraph, ArcMap> RVisitor;
468 468
    RVisitor rvisitor(rgraph, cutMap, cutNum);
469 469

	
470 470
    DfsVisit<RDigraph, RVisitor> rdfs(rgraph, rvisitor);
471 471

	
472 472
    rdfs.init();
473 473
    for (RIterator it = nodes.rbegin(); it != nodes.rend(); ++it) {
474 474
      if (!rdfs.reached(*it)) {
475 475
        rdfs.addSource(*it);
476 476
        rdfs.start();
477 477
      }
478 478
    }
479 479
    return cutNum;
480 480
  }
481 481

	
482 482
  namespace _connectivity_bits {
483 483

	
484 484
    template <typename Digraph>
485 485
    class CountBiNodeConnectedComponentsVisitor : public DfsVisitor<Digraph> {
486 486
    public:
487 487
      typedef typename Digraph::Node Node;
488 488
      typedef typename Digraph::Arc Arc;
489 489
      typedef typename Digraph::Edge Edge;
490 490

	
491 491
      CountBiNodeConnectedComponentsVisitor(const Digraph& graph, int &compNum)
492 492
        : _graph(graph), _compNum(compNum),
493 493
          _numMap(graph), _retMap(graph), _predMap(graph), _num(0) {}
494 494

	
495 495
      void start(const Node& node) {
496 496
        _predMap.set(node, INVALID);
497 497
      }
498 498

	
499 499
      void reach(const Node& node) {
500 500
        _numMap.set(node, _num);
501 501
        _retMap.set(node, _num);
502 502
        ++_num;
503 503
      }
504 504

	
505 505
      void discover(const Arc& edge) {
506 506
        _predMap.set(_graph.target(edge), _graph.source(edge));
507 507
      }
508 508

	
509 509
      void examine(const Arc& edge) {
510 510
        if (_graph.source(edge) == _graph.target(edge) &&
511 511
            _graph.direction(edge)) {
512 512
          ++_compNum;
513 513
          return;
514 514
        }
515 515
        if (_predMap[_graph.source(edge)] == _graph.target(edge)) {
516 516
          return;
517 517
        }
518 518
        if (_retMap[_graph.source(edge)] > _numMap[_graph.target(edge)]) {
519 519
          _retMap.set(_graph.source(edge), _numMap[_graph.target(edge)]);
520 520
        }
521 521
      }
522 522

	
523 523
      void backtrack(const Arc& edge) {
524 524
        if (_retMap[_graph.source(edge)] > _retMap[_graph.target(edge)]) {
525 525
          _retMap.set(_graph.source(edge), _retMap[_graph.target(edge)]);
526 526
        }
527 527
        if (_numMap[_graph.source(edge)] <= _retMap[_graph.target(edge)]) {
528 528
          ++_compNum;
529 529
        }
530 530
      }
531 531

	
532 532
    private:
533 533
      const Digraph& _graph;
534 534
      int& _compNum;
535 535

	
536 536
      typename Digraph::template NodeMap<int> _numMap;
537 537
      typename Digraph::template NodeMap<int> _retMap;
538 538
      typename Digraph::template NodeMap<Node> _predMap;
539 539
      int _num;
540 540
    };
541 541

	
542 542
    template <typename Digraph, typename ArcMap>
543 543
    class BiNodeConnectedComponentsVisitor : public DfsVisitor<Digraph> {
544 544
    public:
545 545
      typedef typename Digraph::Node Node;
546 546
      typedef typename Digraph::Arc Arc;
547 547
      typedef typename Digraph::Edge Edge;
548 548

	
549 549
      BiNodeConnectedComponentsVisitor(const Digraph& graph,
550 550
                                       ArcMap& compMap, int &compNum)
551 551
        : _graph(graph), _compMap(compMap), _compNum(compNum),
552 552
          _numMap(graph), _retMap(graph), _predMap(graph), _num(0) {}
553 553

	
554 554
      void start(const Node& node) {
555 555
        _predMap.set(node, INVALID);
556 556
      }
557 557

	
558 558
      void reach(const Node& node) {
559 559
        _numMap.set(node, _num);
560 560
        _retMap.set(node, _num);
561 561
        ++_num;
562 562
      }
563 563

	
564 564
      void discover(const Arc& edge) {
565 565
        Node target = _graph.target(edge);
566 566
        _predMap.set(target, edge);
567 567
        _edgeStack.push(edge);
568 568
      }
569 569

	
570 570
      void examine(const Arc& edge) {
571 571
        Node source = _graph.source(edge);
572 572
        Node target = _graph.target(edge);
573 573
        if (source == target && _graph.direction(edge)) {
574 574
          _compMap.set(edge, _compNum);
575 575
          ++_compNum;
576 576
          return;
577 577
        }
578 578
        if (_numMap[target] < _numMap[source]) {
579 579
          if (_predMap[source] != _graph.oppositeArc(edge)) {
580 580
            _edgeStack.push(edge);
581 581
          }
582 582
        }
583 583
        if (_predMap[source] != INVALID &&
584 584
            target == _graph.source(_predMap[source])) {
585 585
          return;
586 586
        }
587 587
        if (_retMap[source] > _numMap[target]) {
588 588
          _retMap.set(source, _numMap[target]);
589 589
        }
590 590
      }
591 591

	
592 592
      void backtrack(const Arc& edge) {
593 593
        Node source = _graph.source(edge);
594 594
        Node target = _graph.target(edge);
595 595
        if (_retMap[source] > _retMap[target]) {
596 596
          _retMap.set(source, _retMap[target]);
597 597
        }
598 598
        if (_numMap[source] <= _retMap[target]) {
599 599
          while (_edgeStack.top() != edge) {
600 600
            _compMap.set(_edgeStack.top(), _compNum);
601 601
            _edgeStack.pop();
602 602
          }
603 603
          _compMap.set(edge, _compNum);
604 604
          _edgeStack.pop();
605 605
          ++_compNum;
606 606
        }
607 607
      }
608 608

	
609 609
    private:
610 610
      const Digraph& _graph;
611 611
      ArcMap& _compMap;
612 612
      int& _compNum;
613 613

	
614 614
      typename Digraph::template NodeMap<int> _numMap;
615 615
      typename Digraph::template NodeMap<int> _retMap;
616 616
      typename Digraph::template NodeMap<Arc> _predMap;
617 617
      std::stack<Edge> _edgeStack;
618 618
      int _num;
619 619
    };
620 620

	
621 621

	
622 622
    template <typename Digraph, typename NodeMap>
623 623
    class BiNodeConnectedCutNodesVisitor : public DfsVisitor<Digraph> {
624 624
    public:
625 625
      typedef typename Digraph::Node Node;
626 626
      typedef typename Digraph::Arc Arc;
627 627
      typedef typename Digraph::Edge Edge;
628 628

	
629 629
      BiNodeConnectedCutNodesVisitor(const Digraph& graph, NodeMap& cutMap,
630 630
                                     int& cutNum)
631 631
        : _graph(graph), _cutMap(cutMap), _cutNum(cutNum),
632 632
          _numMap(graph), _retMap(graph), _predMap(graph), _num(0) {}
633 633

	
634 634
      void start(const Node& node) {
635 635
        _predMap.set(node, INVALID);
636 636
        rootCut = false;
637 637
      }
638 638

	
639 639
      void reach(const Node& node) {
640 640
        _numMap.set(node, _num);
641 641
        _retMap.set(node, _num);
642 642
        ++_num;
643 643
      }
644 644

	
645 645
      void discover(const Arc& edge) {
646 646
        _predMap.set(_graph.target(edge), _graph.source(edge));
647 647
      }
648 648

	
649 649
      void examine(const Arc& edge) {
650 650
        if (_graph.source(edge) == _graph.target(edge) &&
651 651
            _graph.direction(edge)) {
652 652
          if (!_cutMap[_graph.source(edge)]) {
653 653
            _cutMap.set(_graph.source(edge), true);
654 654
            ++_cutNum;
655 655
          }
656 656
          return;
657 657
        }
658 658
        if (_predMap[_graph.source(edge)] == _graph.target(edge)) return;
659 659
        if (_retMap[_graph.source(edge)] > _numMap[_graph.target(edge)]) {
660 660
          _retMap.set(_graph.source(edge), _numMap[_graph.target(edge)]);
661 661
        }
662 662
      }
663 663

	
664 664
      void backtrack(const Arc& edge) {
665 665
        if (_retMap[_graph.source(edge)] > _retMap[_graph.target(edge)]) {
666 666
          _retMap.set(_graph.source(edge), _retMap[_graph.target(edge)]);
667 667
        }
668 668
        if (_numMap[_graph.source(edge)] <= _retMap[_graph.target(edge)]) {
669 669
          if (_predMap[_graph.source(edge)] != INVALID) {
670 670
            if (!_cutMap[_graph.source(edge)]) {
671 671
              _cutMap.set(_graph.source(edge), true);
672 672
              ++_cutNum;
673 673
            }
674 674
          } else if (rootCut) {
675 675
            if (!_cutMap[_graph.source(edge)]) {
676 676
              _cutMap.set(_graph.source(edge), true);
677 677
              ++_cutNum;
678 678
            }
679 679
          } else {
680 680
            rootCut = true;
681 681
          }
682 682
        }
683 683
      }
684 684

	
685 685
    private:
686 686
      const Digraph& _graph;
687 687
      NodeMap& _cutMap;
688 688
      int& _cutNum;
689 689

	
690 690
      typename Digraph::template NodeMap<int> _numMap;
691 691
      typename Digraph::template NodeMap<int> _retMap;
692 692
      typename Digraph::template NodeMap<Node> _predMap;
693 693
      std::stack<Edge> _edgeStack;
694 694
      int _num;
695 695
      bool rootCut;
696 696
    };
697 697

	
698 698
  }
699 699

	
700 700
  template <typename Graph>
701 701
  int countBiNodeConnectedComponents(const Graph& graph);
702 702

	
703 703
  /// \ingroup connectivity
704 704
  ///
705 705
  /// \brief Checks the graph is bi-node-connected.
706 706
  ///
707 707
  /// This function checks that the undirected graph is bi-node-connected
708 708
  /// graph. The graph is bi-node-connected if any two undirected edge is
709 709
  /// on same circle.
710 710
  ///
711 711
  /// \param graph The graph.
712 712
  /// \return %True when the graph bi-node-connected.
713 713
  template <typename Graph>
714 714
  bool biNodeConnected(const Graph& graph) {
715 715
    return countBiNodeConnectedComponents(graph) <= 1;
716 716
  }
717 717

	
718 718
  /// \ingroup connectivity
719 719
  ///
720 720
  /// \brief Count the biconnected components.
721 721
  ///
722 722
  /// This function finds the bi-node-connected components in an undirected
723 723
  /// graph. The biconnected components are the classes of an equivalence
724 724
  /// relation on the undirected edges. Two undirected edge is in relationship
725 725
  /// when they are on same circle.
726 726
  ///
727 727
  /// \param graph The graph.
728 728
  /// \return The number of components.
729 729
  template <typename Graph>
730 730
  int countBiNodeConnectedComponents(const Graph& graph) {
731 731
    checkConcept<concepts::Graph, Graph>();
732 732
    typedef typename Graph::NodeIt NodeIt;
733 733

	
734 734
    using namespace _connectivity_bits;
735 735

	
736 736
    typedef CountBiNodeConnectedComponentsVisitor<Graph> Visitor;
737 737

	
738 738
    int compNum = 0;
739 739
    Visitor visitor(graph, compNum);
740 740

	
741 741
    DfsVisit<Graph, Visitor> dfs(graph, visitor);
742 742
    dfs.init();
743 743

	
744 744
    for (NodeIt it(graph); it != INVALID; ++it) {
745 745
      if (!dfs.reached(it)) {
746 746
        dfs.addSource(it);
747 747
        dfs.start();
748 748
      }
749 749
    }
750 750
    return compNum;
751 751
  }
752 752

	
753 753
  /// \ingroup connectivity
754 754
  ///
755 755
  /// \brief Find the bi-node-connected components.
756 756
  ///
757 757
  /// This function finds the bi-node-connected components in an undirected
758 758
  /// graph. The bi-node-connected components are the classes of an equivalence
759 759
  /// relation on the undirected edges. Two undirected edge are in relationship
760 760
  /// when they are on same circle.
761 761
  ///
762 762
  /// \param graph The graph.
763 763
  /// \retval compMap A writable uedge map. The values will be set from 0
764 764
  /// to the number of the biconnected components minus one. Each values
765 765
  /// of the map will be set exactly once, the values of a certain component
766 766
  /// will be set continuously.
767 767
  /// \return The number of components.
768 768
  ///
769 769
  template <typename Graph, typename EdgeMap>
770 770
  int biNodeConnectedComponents(const Graph& graph,
771 771
                                EdgeMap& compMap) {
772 772
    checkConcept<concepts::Graph, Graph>();
773 773
    typedef typename Graph::NodeIt NodeIt;
774 774
    typedef typename Graph::Edge Edge;
775 775
    checkConcept<concepts::WriteMap<Edge, int>, EdgeMap>();
776 776

	
777 777
    using namespace _connectivity_bits;
778 778

	
779 779
    typedef BiNodeConnectedComponentsVisitor<Graph, EdgeMap> Visitor;
780 780

	
781 781
    int compNum = 0;
782 782
    Visitor visitor(graph, compMap, compNum);
783 783

	
784 784
    DfsVisit<Graph, Visitor> dfs(graph, visitor);
785 785
    dfs.init();
786 786

	
787 787
    for (NodeIt it(graph); it != INVALID; ++it) {
788 788
      if (!dfs.reached(it)) {
789 789
        dfs.addSource(it);
790 790
        dfs.start();
791 791
      }
792 792
    }
793 793
    return compNum;
794 794
  }
795 795

	
796 796
  /// \ingroup connectivity
797 797
  ///
798 798
  /// \brief Find the bi-node-connected cut nodes.
799 799
  ///
800 800
  /// This function finds the bi-node-connected cut nodes in an undirected
801 801
  /// graph. The bi-node-connected components are the classes of an equivalence
802 802
  /// relation on the undirected edges. Two undirected edges are in
803 803
  /// relationship when they are on same circle. The biconnected components
804 804
  /// are separted by nodes which are the cut nodes of the components.
805 805
  ///
806 806
  /// \param graph The graph.
807 807
  /// \retval cutMap A writable edge map. The values will be set true when
808 808
  /// the node separate two or more components.
809 809
  /// \return The number of the cut nodes.
810 810
  template <typename Graph, typename NodeMap>
811 811
  int biNodeConnectedCutNodes(const Graph& graph, NodeMap& cutMap) {
812 812
    checkConcept<concepts::Graph, Graph>();
813 813
    typedef typename Graph::Node Node;
814 814
    typedef typename Graph::NodeIt NodeIt;
815 815
    checkConcept<concepts::WriteMap<Node, bool>, NodeMap>();
816 816

	
817 817
    using namespace _connectivity_bits;
818 818

	
819 819
    typedef BiNodeConnectedCutNodesVisitor<Graph, NodeMap> Visitor;
820 820

	
821 821
    int cutNum = 0;
822 822
    Visitor visitor(graph, cutMap, cutNum);
823 823

	
824 824
    DfsVisit<Graph, Visitor> dfs(graph, visitor);
825 825
    dfs.init();
826 826

	
827 827
    for (NodeIt it(graph); it != INVALID; ++it) {
828 828
      if (!dfs.reached(it)) {
829 829
        dfs.addSource(it);
830 830
        dfs.start();
831 831
      }
832 832
    }
833 833
    return cutNum;
834 834
  }
835 835

	
836 836
  namespace _connectivity_bits {
837 837

	
838 838
    template <typename Digraph>
839 839
    class CountBiEdgeConnectedComponentsVisitor : public DfsVisitor<Digraph> {
840 840
    public:
841 841
      typedef typename Digraph::Node Node;
842 842
      typedef typename Digraph::Arc Arc;
843 843
      typedef typename Digraph::Edge Edge;
844 844

	
845 845
      CountBiEdgeConnectedComponentsVisitor(const Digraph& graph, int &compNum)
846 846
        : _graph(graph), _compNum(compNum),
847 847
          _numMap(graph), _retMap(graph), _predMap(graph), _num(0) {}
848 848

	
849 849
      void start(const Node& node) {
850 850
        _predMap.set(node, INVALID);
851 851
      }
852 852

	
853 853
      void reach(const Node& node) {
854 854
        _numMap.set(node, _num);
855 855
        _retMap.set(node, _num);
856 856
        ++_num;
857 857
      }
858 858

	
859 859
      void leave(const Node& node) {
860 860
        if (_numMap[node] <= _retMap[node]) {
861 861
          ++_compNum;
862 862
        }
863 863
      }
864 864

	
865 865
      void discover(const Arc& edge) {
866 866
        _predMap.set(_graph.target(edge), edge);
867 867
      }
868 868

	
869 869
      void examine(const Arc& edge) {
870 870
        if (_predMap[_graph.source(edge)] == _graph.oppositeArc(edge)) {
871 871
          return;
872 872
        }
873 873
        if (_retMap[_graph.source(edge)] > _retMap[_graph.target(edge)]) {
874 874
          _retMap.set(_graph.source(edge), _retMap[_graph.target(edge)]);
875 875
        }
876 876
      }
877 877

	
878 878
      void backtrack(const Arc& edge) {
879 879
        if (_retMap[_graph.source(edge)] > _retMap[_graph.target(edge)]) {
880 880
          _retMap.set(_graph.source(edge), _retMap[_graph.target(edge)]);
881 881
        }
882 882
      }
883 883

	
884 884
    private:
885 885
      const Digraph& _graph;
886 886
      int& _compNum;
887 887

	
888 888
      typename Digraph::template NodeMap<int> _numMap;
889 889
      typename Digraph::template NodeMap<int> _retMap;
890 890
      typename Digraph::template NodeMap<Arc> _predMap;
891 891
      int _num;
892 892
    };
893 893

	
894 894
    template <typename Digraph, typename NodeMap>
895 895
    class BiEdgeConnectedComponentsVisitor : public DfsVisitor<Digraph> {
896 896
    public:
897 897
      typedef typename Digraph::Node Node;
898 898
      typedef typename Digraph::Arc Arc;
899 899
      typedef typename Digraph::Edge Edge;
900 900

	
901 901
      BiEdgeConnectedComponentsVisitor(const Digraph& graph,
902 902
                                       NodeMap& compMap, int &compNum)
903 903
        : _graph(graph), _compMap(compMap), _compNum(compNum),
904 904
          _numMap(graph), _retMap(graph), _predMap(graph), _num(0) {}
905 905

	
906 906
      void start(const Node& node) {
907 907
        _predMap.set(node, INVALID);
908 908
      }
909 909

	
910 910
      void reach(const Node& node) {
911 911
        _numMap.set(node, _num);
912 912
        _retMap.set(node, _num);
913 913
        _nodeStack.push(node);
914 914
        ++_num;
915 915
      }
916 916

	
917 917
      void leave(const Node& node) {
918 918
        if (_numMap[node] <= _retMap[node]) {
919 919
          while (_nodeStack.top() != node) {
920 920
            _compMap.set(_nodeStack.top(), _compNum);
921 921
            _nodeStack.pop();
922 922
          }
923 923
          _compMap.set(node, _compNum);
924 924
          _nodeStack.pop();
925 925
          ++_compNum;
926 926
        }
927 927
      }
928 928

	
929 929
      void discover(const Arc& edge) {
930 930
        _predMap.set(_graph.target(edge), edge);
931 931
      }
932 932

	
933 933
      void examine(const Arc& edge) {
934 934
        if (_predMap[_graph.source(edge)] == _graph.oppositeArc(edge)) {
935 935
          return;
936 936
        }
937 937
        if (_retMap[_graph.source(edge)] > _retMap[_graph.target(edge)]) {
938 938
          _retMap.set(_graph.source(edge), _retMap[_graph.target(edge)]);
939 939
        }
940 940
      }
941 941

	
942 942
      void backtrack(const Arc& edge) {
943 943
        if (_retMap[_graph.source(edge)] > _retMap[_graph.target(edge)]) {
944 944
          _retMap.set(_graph.source(edge), _retMap[_graph.target(edge)]);
945 945
        }
946 946
      }
947 947

	
948 948
    private:
949 949
      const Digraph& _graph;
950 950
      NodeMap& _compMap;
951 951
      int& _compNum;
952 952

	
953 953
      typename Digraph::template NodeMap<int> _numMap;
954 954
      typename Digraph::template NodeMap<int> _retMap;
955 955
      typename Digraph::template NodeMap<Arc> _predMap;
956 956
      std::stack<Node> _nodeStack;
957 957
      int _num;
958 958
    };
959 959

	
960 960

	
961 961
    template <typename Digraph, typename ArcMap>
962 962
    class BiEdgeConnectedCutEdgesVisitor : public DfsVisitor<Digraph> {
963 963
    public:
964 964
      typedef typename Digraph::Node Node;
965 965
      typedef typename Digraph::Arc Arc;
966 966
      typedef typename Digraph::Edge Edge;
967 967

	
968 968
      BiEdgeConnectedCutEdgesVisitor(const Digraph& graph,
969 969
                                     ArcMap& cutMap, int &cutNum)
970 970
        : _graph(graph), _cutMap(cutMap), _cutNum(cutNum),
971 971
          _numMap(graph), _retMap(graph), _predMap(graph), _num(0) {}
972 972

	
973 973
      void start(const Node& node) {
974 974
        _predMap[node] = INVALID;
975 975
      }
976 976

	
977 977
      void reach(const Node& node) {
978 978
        _numMap.set(node, _num);
979 979
        _retMap.set(node, _num);
980 980
        ++_num;
981 981
      }
982 982

	
983 983
      void leave(const Node& node) {
984 984
        if (_numMap[node] <= _retMap[node]) {
985 985
          if (_predMap[node] != INVALID) {
986 986
            _cutMap.set(_predMap[node], true);
987 987
            ++_cutNum;
988 988
          }
989 989
        }
990 990
      }
991 991

	
992 992
      void discover(const Arc& edge) {
993 993
        _predMap.set(_graph.target(edge), edge);
994 994
      }
995 995

	
996 996
      void examine(const Arc& edge) {
997 997
        if (_predMap[_graph.source(edge)] == _graph.oppositeArc(edge)) {
998 998
          return;
999 999
        }
1000 1000
        if (_retMap[_graph.source(edge)] > _retMap[_graph.target(edge)]) {
1001 1001
          _retMap.set(_graph.source(edge), _retMap[_graph.target(edge)]);
1002 1002
        }
1003 1003
      }
1004 1004

	
1005 1005
      void backtrack(const Arc& edge) {
1006 1006
        if (_retMap[_graph.source(edge)] > _retMap[_graph.target(edge)]) {
1007 1007
          _retMap.set(_graph.source(edge), _retMap[_graph.target(edge)]);
1008 1008
        }
1009 1009
      }
1010 1010

	
1011 1011
    private:
1012 1012
      const Digraph& _graph;
1013 1013
      ArcMap& _cutMap;
1014 1014
      int& _cutNum;
1015 1015

	
1016 1016
      typename Digraph::template NodeMap<int> _numMap;
1017 1017
      typename Digraph::template NodeMap<int> _retMap;
1018 1018
      typename Digraph::template NodeMap<Arc> _predMap;
1019 1019
      int _num;
1020 1020
    };
1021 1021
  }
1022 1022

	
1023 1023
  template <typename Graph>
1024 1024
  int countBiEdgeConnectedComponents(const Graph& graph);
1025 1025

	
1026 1026
  /// \ingroup connectivity
1027 1027
  ///
1028 1028
  /// \brief Checks that the graph is bi-edge-connected.
1029 1029
  ///
1030 1030
  /// This function checks that the graph is bi-edge-connected. The undirected
1031 1031
  /// graph is bi-edge-connected when any two nodes are connected with two
1032 1032
  /// edge-disjoint paths.
1033 1033
  ///
1034 1034
  /// \param graph The undirected graph.
1035 1035
  /// \return The number of components.
1036 1036
  template <typename Graph>
1037 1037
  bool biEdgeConnected(const Graph& graph) {
1038 1038
    return countBiEdgeConnectedComponents(graph) <= 1;
1039 1039
  }
1040 1040

	
1041 1041
  /// \ingroup connectivity
1042 1042
  ///
1043 1043
  /// \brief Count the bi-edge-connected components.
1044 1044
  ///
1045 1045
  /// This function count the bi-edge-connected components in an undirected
1046 1046
  /// graph. The bi-edge-connected components are the classes of an equivalence
1047 1047
  /// relation on the nodes. Two nodes are in relationship when they are
1048 1048
  /// connected with at least two edge-disjoint paths.
1049 1049
  ///
1050 1050
  /// \param graph The undirected graph.
1051 1051
  /// \return The number of components.
1052 1052
  template <typename Graph>
1053 1053
  int countBiEdgeConnectedComponents(const Graph& graph) {
1054 1054
    checkConcept<concepts::Graph, Graph>();
1055 1055
    typedef typename Graph::NodeIt NodeIt;
1056 1056

	
1057 1057
    using namespace _connectivity_bits;
1058 1058

	
1059 1059
    typedef CountBiEdgeConnectedComponentsVisitor<Graph> Visitor;
1060 1060

	
1061 1061
    int compNum = 0;
1062 1062
    Visitor visitor(graph, compNum);
1063 1063

	
1064 1064
    DfsVisit<Graph, Visitor> dfs(graph, visitor);
1065 1065
    dfs.init();
1066 1066

	
1067 1067
    for (NodeIt it(graph); it != INVALID; ++it) {
1068 1068
      if (!dfs.reached(it)) {
1069 1069
        dfs.addSource(it);
1070 1070
        dfs.start();
1071 1071
      }
1072 1072
    }
1073 1073
    return compNum;
1074 1074
  }
1075 1075

	
1076 1076
  /// \ingroup connectivity
1077 1077
  ///
1078 1078
  /// \brief Find the bi-edge-connected components.
1079 1079
  ///
1080 1080
  /// This function finds the bi-edge-connected components in an undirected
1081 1081
  /// graph. The bi-edge-connected components are the classes of an equivalence
1082 1082
  /// relation on the nodes. Two nodes are in relationship when they are
1083 1083
  /// connected at least two edge-disjoint paths.
1084 1084
  ///
1085 1085
  /// \param graph The graph.
1086 1086
  /// \retval compMap A writable node map. The values will be set from 0 to
1087 1087
  /// the number of the biconnected components minus one. Each values
1088 1088
  /// of the map will be set exactly once, the values of a certain component
1089 1089
  /// will be set continuously.
1090 1090
  /// \return The number of components.
1091 1091
  ///
1092 1092
  template <typename Graph, typename NodeMap>
1093 1093
  int biEdgeConnectedComponents(const Graph& graph, NodeMap& compMap) {
1094 1094
    checkConcept<concepts::Graph, Graph>();
1095 1095
    typedef typename Graph::NodeIt NodeIt;
1096 1096
    typedef typename Graph::Node Node;
1097 1097
    checkConcept<concepts::WriteMap<Node, int>, NodeMap>();
1098 1098

	
1099 1099
    using namespace _connectivity_bits;
1100 1100

	
1101 1101
    typedef BiEdgeConnectedComponentsVisitor<Graph, NodeMap> Visitor;
1102 1102

	
1103 1103
    int compNum = 0;
1104 1104
    Visitor visitor(graph, compMap, compNum);
1105 1105

	
1106 1106
    DfsVisit<Graph, Visitor> dfs(graph, visitor);
1107 1107
    dfs.init();
1108 1108

	
1109 1109
    for (NodeIt it(graph); it != INVALID; ++it) {
1110 1110
      if (!dfs.reached(it)) {
1111 1111
        dfs.addSource(it);
1112 1112
        dfs.start();
1113 1113
      }
1114 1114
    }
1115 1115
    return compNum;
1116 1116
  }
1117 1117

	
1118 1118
  /// \ingroup connectivity
1119 1119
  ///
1120 1120
  /// \brief Find the bi-edge-connected cut edges.
1121 1121
  ///
1122 1122
  /// This function finds the bi-edge-connected components in an undirected
1123 1123
  /// graph. The bi-edge-connected components are the classes of an equivalence
1124 1124
  /// relation on the nodes. Two nodes are in relationship when they are
1125 1125
  /// connected with at least two edge-disjoint paths. The bi-edge-connected
1126 1126
  /// components are separted by edges which are the cut edges of the
1127 1127
  /// components.
1128 1128
  ///
1129 1129
  /// \param graph The graph.
1130 1130
  /// \retval cutMap A writable node map. The values will be set true when the
1131 1131
  /// edge is a cut edge.
1132 1132
  /// \return The number of cut edges.
1133 1133
  template <typename Graph, typename EdgeMap>
1134 1134
  int biEdgeConnectedCutEdges(const Graph& graph, EdgeMap& cutMap) {
1135 1135
    checkConcept<concepts::Graph, Graph>();
1136 1136
    typedef typename Graph::NodeIt NodeIt;
1137 1137
    typedef typename Graph::Edge Edge;
1138 1138
    checkConcept<concepts::WriteMap<Edge, bool>, EdgeMap>();
1139 1139

	
1140 1140
    using namespace _connectivity_bits;
1141 1141

	
1142 1142
    typedef BiEdgeConnectedCutEdgesVisitor<Graph, EdgeMap> Visitor;
1143 1143

	
1144 1144
    int cutNum = 0;
1145 1145
    Visitor visitor(graph, cutMap, cutNum);
1146 1146

	
1147 1147
    DfsVisit<Graph, Visitor> dfs(graph, visitor);
1148 1148
    dfs.init();
1149 1149

	
1150 1150
    for (NodeIt it(graph); it != INVALID; ++it) {
1151 1151
      if (!dfs.reached(it)) {
1152 1152
        dfs.addSource(it);
1153 1153
        dfs.start();
1154 1154
      }
1155 1155
    }
1156 1156
    return cutNum;
1157 1157
  }
1158 1158

	
1159 1159

	
1160 1160
  namespace _connectivity_bits {
1161 1161

	
1162 1162
    template <typename Digraph, typename IntNodeMap>
1163 1163
    class TopologicalSortVisitor : public DfsVisitor<Digraph> {
1164 1164
    public:
1165 1165
      typedef typename Digraph::Node Node;
1166 1166
      typedef typename Digraph::Arc edge;
1167 1167

	
1168 1168
      TopologicalSortVisitor(IntNodeMap& order, int num)
1169 1169
        : _order(order), _num(num) {}
1170 1170

	
1171 1171
      void leave(const Node& node) {
1172 1172
        _order.set(node, --_num);
1173 1173
      }
1174 1174

	
1175 1175
    private:
1176 1176
      IntNodeMap& _order;
1177 1177
      int _num;
1178 1178
    };
1179 1179

	
1180 1180
  }
1181 1181

	
1182 1182
  /// \ingroup connectivity
1183 1183
  ///
1184 1184
  /// \brief Sort the nodes of a DAG into topolgical order.
1185 1185
  ///
1186 1186
  /// Sort the nodes of a DAG into topolgical order.
1187 1187
  ///
1188 1188
  /// \param graph The graph. It must be directed and acyclic.
1189 1189
  /// \retval order A writable node map. The values will be set from 0 to
1190 1190
  /// the number of the nodes in the graph minus one. Each values of the map
1191 1191
  /// will be set exactly once, the values  will be set descending order.
1192 1192
  ///
1193 1193
  /// \see checkedTopologicalSort
1194 1194
  /// \see dag
1195 1195
  template <typename Digraph, typename NodeMap>
1196 1196
  void topologicalSort(const Digraph& graph, NodeMap& order) {
1197 1197
    using namespace _connectivity_bits;
1198 1198

	
1199 1199
    checkConcept<concepts::Digraph, Digraph>();
1200 1200
    checkConcept<concepts::WriteMap<typename Digraph::Node, int>, NodeMap>();
1201 1201

	
1202 1202
    typedef typename Digraph::Node Node;
1203 1203
    typedef typename Digraph::NodeIt NodeIt;
1204 1204
    typedef typename Digraph::Arc Arc;
1205 1205

	
1206 1206
    TopologicalSortVisitor<Digraph, NodeMap>
1207 1207
      visitor(order, countNodes(graph));
1208 1208

	
1209 1209
    DfsVisit<Digraph, TopologicalSortVisitor<Digraph, NodeMap> >
1210 1210
      dfs(graph, visitor);
1211 1211

	
1212 1212
    dfs.init();
1213 1213
    for (NodeIt it(graph); it != INVALID; ++it) {
1214 1214
      if (!dfs.reached(it)) {
1215 1215
        dfs.addSource(it);
1216 1216
        dfs.start();
1217 1217
      }
1218 1218
    }
1219 1219
  }
1220 1220

	
1221 1221
  /// \ingroup connectivity
1222 1222
  ///
1223 1223
  /// \brief Sort the nodes of a DAG into topolgical order.
1224 1224
  ///
1225 1225
  /// Sort the nodes of a DAG into topolgical order. It also checks
1226 1226
  /// that the given graph is DAG.
1227 1227
  ///
1228
  /// \param graph The graph. It must be directed and acyclic.
1228
  /// \param digraph The graph. It must be directed and acyclic.
1229 1229
  /// \retval order A readable - writable node map. The values will be set
1230 1230
  /// from 0 to the number of the nodes in the graph minus one. Each values
1231 1231
  /// of the map will be set exactly once, the values will be set descending
1232 1232
  /// order.
1233 1233
  /// \return %False when the graph is not DAG.
1234 1234
  ///
1235 1235
  /// \see topologicalSort
1236 1236
  /// \see dag
1237 1237
  template <typename Digraph, typename NodeMap>
1238 1238
  bool checkedTopologicalSort(const Digraph& digraph, NodeMap& order) {
1239 1239
    using namespace _connectivity_bits;
1240 1240

	
1241 1241
    checkConcept<concepts::Digraph, Digraph>();
1242 1242
    checkConcept<concepts::ReadWriteMap<typename Digraph::Node, int>,
1243 1243
      NodeMap>();
1244 1244

	
1245 1245
    typedef typename Digraph::Node Node;
1246 1246
    typedef typename Digraph::NodeIt NodeIt;
1247 1247
    typedef typename Digraph::Arc Arc;
1248 1248

	
1249 1249
    for (NodeIt it(digraph); it != INVALID; ++it) {
1250 1250
      order.set(it, -1);
1251 1251
    }
1252 1252

	
1253 1253
    TopologicalSortVisitor<Digraph, NodeMap>
1254 1254
      visitor(order, countNodes(digraph));
1255 1255

	
1256 1256
    DfsVisit<Digraph, TopologicalSortVisitor<Digraph, NodeMap> >
1257 1257
      dfs(digraph, visitor);
1258 1258

	
1259 1259
    dfs.init();
1260 1260
    for (NodeIt it(digraph); it != INVALID; ++it) {
1261 1261
      if (!dfs.reached(it)) {
1262 1262
        dfs.addSource(it);
1263 1263
        while (!dfs.emptyQueue()) {
1264 1264
           Arc arc = dfs.nextArc();
1265 1265
           Node target = digraph.target(arc);
1266 1266
           if (dfs.reached(target) && order[target] == -1) {
1267 1267
             return false;
1268 1268
           }
1269 1269
           dfs.processNextArc();
1270 1270
         }
1271 1271
      }
1272 1272
    }
1273 1273
    return true;
1274 1274
  }
1275 1275

	
1276 1276
  /// \ingroup connectivity
1277 1277
  ///
1278 1278
  /// \brief Check that the given directed graph is a DAG.
1279 1279
  ///
1280 1280
  /// Check that the given directed graph is a DAG. The DAG is
1281 1281
  /// an Directed Acyclic Digraph.
1282 1282
  /// \return %False when the graph is not DAG.
1283 1283
  /// \see acyclic
1284 1284
  template <typename Digraph>
1285 1285
  bool dag(const Digraph& digraph) {
1286 1286

	
1287 1287
    checkConcept<concepts::Digraph, Digraph>();
1288 1288

	
1289 1289
    typedef typename Digraph::Node Node;
1290 1290
    typedef typename Digraph::NodeIt NodeIt;
1291 1291
    typedef typename Digraph::Arc Arc;
1292 1292

	
1293 1293
    typedef typename Digraph::template NodeMap<bool> ProcessedMap;
1294 1294

	
1295 1295
    typename Dfs<Digraph>::template SetProcessedMap<ProcessedMap>::
1296 1296
      Create dfs(digraph);
1297 1297

	
1298 1298
    ProcessedMap processed(digraph);
1299 1299
    dfs.processedMap(processed);
1300 1300

	
1301 1301
    dfs.init();
1302 1302
    for (NodeIt it(digraph); it != INVALID; ++it) {
1303 1303
      if (!dfs.reached(it)) {
1304 1304
        dfs.addSource(it);
1305 1305
        while (!dfs.emptyQueue()) {
1306 1306
          Arc edge = dfs.nextArc();
1307 1307
          Node target = digraph.target(edge);
1308 1308
          if (dfs.reached(target) && !processed[target]) {
1309 1309
            return false;
1310 1310
          }
1311 1311
          dfs.processNextArc();
1312 1312
        }
1313 1313
      }
1314 1314
    }
1315 1315
    return true;
1316 1316
  }
1317 1317

	
1318 1318
  /// \ingroup connectivity
1319 1319
  ///
1320 1320
  /// \brief Check that the given undirected graph is acyclic.
1321 1321
  ///
1322 1322
  /// Check that the given undirected graph acyclic.
1323 1323
  /// \param graph The undirected graph.
1324 1324
  /// \return %True when there is no circle in the graph.
1325 1325
  /// \see dag
1326 1326
  template <typename Graph>
1327 1327
  bool acyclic(const Graph& graph) {
1328 1328
    checkConcept<concepts::Graph, Graph>();
1329 1329
    typedef typename Graph::Node Node;
1330 1330
    typedef typename Graph::NodeIt NodeIt;
1331 1331
    typedef typename Graph::Arc Arc;
1332 1332
    Dfs<Graph> dfs(graph);
1333 1333
    dfs.init();
1334 1334
    for (NodeIt it(graph); it != INVALID; ++it) {
1335 1335
      if (!dfs.reached(it)) {
1336 1336
        dfs.addSource(it);
1337 1337
        while (!dfs.emptyQueue()) {
1338 1338
          Arc edge = dfs.nextArc();
1339 1339
          Node source = graph.source(edge);
1340 1340
          Node target = graph.target(edge);
1341 1341
          if (dfs.reached(target) &&
1342 1342
              dfs.predArc(source) != graph.oppositeArc(edge)) {
1343 1343
            return false;
1344 1344
          }
1345 1345
          dfs.processNextArc();
1346 1346
        }
1347 1347
      }
1348 1348
    }
1349 1349
    return true;
1350 1350
  }
1351 1351

	
1352 1352
  /// \ingroup connectivity
1353 1353
  ///
1354 1354
  /// \brief Check that the given undirected graph is tree.
1355 1355
  ///
1356 1356
  /// Check that the given undirected graph is tree.
1357 1357
  /// \param graph The undirected graph.
1358 1358
  /// \return %True when the graph is acyclic and connected.
1359 1359
  template <typename Graph>
1360 1360
  bool tree(const Graph& graph) {
1361 1361
    checkConcept<concepts::Graph, Graph>();
1362 1362
    typedef typename Graph::Node Node;
1363 1363
    typedef typename Graph::NodeIt NodeIt;
1364 1364
    typedef typename Graph::Arc Arc;
1365 1365
    Dfs<Graph> dfs(graph);
1366 1366
    dfs.init();
1367 1367
    dfs.addSource(NodeIt(graph));
1368 1368
    while (!dfs.emptyQueue()) {
1369 1369
      Arc edge = dfs.nextArc();
1370 1370
      Node source = graph.source(edge);
1371 1371
      Node target = graph.target(edge);
1372 1372
      if (dfs.reached(target) &&
1373 1373
          dfs.predArc(source) != graph.oppositeArc(edge)) {
1374 1374
        return false;
1375 1375
      }
1376 1376
      dfs.processNextArc();
1377 1377
    }
1378 1378
    for (NodeIt it(graph); it != INVALID; ++it) {
1379 1379
      if (!dfs.reached(it)) {
1380 1380
        return false;
1381 1381
      }
1382 1382
    }
1383 1383
    return true;
1384 1384
  }
1385 1385

	
1386 1386
  namespace _connectivity_bits {
1387 1387

	
1388 1388
    template <typename Digraph>
1389 1389
    class BipartiteVisitor : public BfsVisitor<Digraph> {
1390 1390
    public:
1391 1391
      typedef typename Digraph::Arc Arc;
1392 1392
      typedef typename Digraph::Node Node;
1393 1393

	
1394 1394
      BipartiteVisitor(const Digraph& graph, bool& bipartite)
1395 1395
        : _graph(graph), _part(graph), _bipartite(bipartite) {}
1396 1396

	
1397 1397
      void start(const Node& node) {
1398 1398
        _part[node] = true;
1399 1399
      }
1400 1400
      void discover(const Arc& edge) {
1401 1401
        _part.set(_graph.target(edge), !_part[_graph.source(edge)]);
1402 1402
      }
1403 1403
      void examine(const Arc& edge) {
1404 1404
        _bipartite = _bipartite &&
1405 1405
          _part[_graph.target(edge)] != _part[_graph.source(edge)];
1406 1406
      }
1407 1407

	
1408 1408
    private:
1409 1409

	
1410 1410
      const Digraph& _graph;
1411 1411
      typename Digraph::template NodeMap<bool> _part;
1412 1412
      bool& _bipartite;
1413 1413
    };
1414 1414

	
1415 1415
    template <typename Digraph, typename PartMap>
1416 1416
    class BipartitePartitionsVisitor : public BfsVisitor<Digraph> {
1417 1417
    public:
1418 1418
      typedef typename Digraph::Arc Arc;
1419 1419
      typedef typename Digraph::Node Node;
1420 1420

	
1421 1421
      BipartitePartitionsVisitor(const Digraph& graph,
1422 1422
                                 PartMap& part, bool& bipartite)
1423 1423
        : _graph(graph), _part(part), _bipartite(bipartite) {}
1424 1424

	
1425 1425
      void start(const Node& node) {
1426 1426
        _part.set(node, true);
1427 1427
      }
1428 1428
      void discover(const Arc& edge) {
1429 1429
        _part.set(_graph.target(edge), !_part[_graph.source(edge)]);
1430 1430
      }
1431 1431
      void examine(const Arc& edge) {
1432 1432
        _bipartite = _bipartite &&
1433 1433
          _part[_graph.target(edge)] != _part[_graph.source(edge)];
1434 1434
      }
1435 1435

	
1436 1436
    private:
1437 1437

	
1438 1438
      const Digraph& _graph;
1439 1439
      PartMap& _part;
1440 1440
      bool& _bipartite;
1441 1441
    };
1442 1442
  }
1443 1443

	
1444 1444
  /// \ingroup connectivity
1445 1445
  ///
1446 1446
  /// \brief Check if the given undirected graph is bipartite or not
1447 1447
  ///
1448 1448
  /// The function checks if the given undirected \c graph graph is bipartite
1449 1449
  /// or not. The \ref Bfs algorithm is used to calculate the result.
1450 1450
  /// \param graph The undirected graph.
1451 1451
  /// \return %True if \c graph is bipartite, %false otherwise.
1452 1452
  /// \sa bipartitePartitions
1453 1453
  template<typename Graph>
1454 1454
  inline bool bipartite(const Graph &graph){
1455 1455
    using namespace _connectivity_bits;
1456 1456

	
1457 1457
    checkConcept<concepts::Graph, Graph>();
1458 1458

	
1459 1459
    typedef typename Graph::NodeIt NodeIt;
1460 1460
    typedef typename Graph::ArcIt ArcIt;
1461 1461

	
1462 1462
    bool bipartite = true;
1463 1463

	
1464 1464
    BipartiteVisitor<Graph>
1465 1465
      visitor(graph, bipartite);
1466 1466
    BfsVisit<Graph, BipartiteVisitor<Graph> >
1467 1467
      bfs(graph, visitor);
1468 1468
    bfs.init();
1469 1469
    for(NodeIt it(graph); it != INVALID; ++it) {
1470 1470
      if(!bfs.reached(it)){
1471 1471
        bfs.addSource(it);
1472 1472
        while (!bfs.emptyQueue()) {
1473 1473
          bfs.processNextNode();
1474 1474
          if (!bipartite) return false;
1475 1475
        }
1476 1476
      }
1477 1477
    }
1478 1478
    return true;
1479 1479
  }
1480 1480

	
1481 1481
  /// \ingroup connectivity
1482 1482
  ///
1483 1483
  /// \brief Check if the given undirected graph is bipartite or not
1484 1484
  ///
1485 1485
  /// The function checks if the given undirected graph is bipartite
1486 1486
  /// or not. The  \ref  Bfs  algorithm  is   used  to  calculate the result.
1487 1487
  /// During the execution, the \c partMap will be set as the two
1488 1488
  /// partitions of the graph.
1489 1489
  /// \param graph The undirected graph.
1490 1490
  /// \retval partMap A writable bool map of nodes. It will be set as the
1491 1491
  /// two partitions of the graph.
1492 1492
  /// \return %True if \c graph is bipartite, %false otherwise.
1493 1493
  template<typename Graph, typename NodeMap>
1494 1494
  inline bool bipartitePartitions(const Graph &graph, NodeMap &partMap){
1495 1495
    using namespace _connectivity_bits;
1496 1496

	
1497 1497
    checkConcept<concepts::Graph, Graph>();
1498 1498

	
1499 1499
    typedef typename Graph::Node Node;
1500 1500
    typedef typename Graph::NodeIt NodeIt;
1501 1501
    typedef typename Graph::ArcIt ArcIt;
1502 1502

	
1503 1503
    bool bipartite = true;
1504 1504

	
1505 1505
    BipartitePartitionsVisitor<Graph, NodeMap>
1506 1506
      visitor(graph, partMap, bipartite);
1507 1507
    BfsVisit<Graph, BipartitePartitionsVisitor<Graph, NodeMap> >
1508 1508
      bfs(graph, visitor);
1509 1509
    bfs.init();
1510 1510
    for(NodeIt it(graph); it != INVALID; ++it) {
1511 1511
      if(!bfs.reached(it)){
1512 1512
        bfs.addSource(it);
1513 1513
        while (!bfs.emptyQueue()) {
1514 1514
          bfs.processNextNode();
1515 1515
          if (!bipartite) return false;
1516 1516
        }
1517 1517
      }
1518 1518
    }
1519 1519
    return true;
1520 1520
  }
1521 1521

	
1522 1522
  /// \brief Returns true when there are not loop edges in the graph.
1523 1523
  ///
1524 1524
  /// Returns true when there are not loop edges in the graph.
1525 1525
  template <typename Digraph>
1526 1526
  bool loopFree(const Digraph& digraph) {
1527 1527
    for (typename Digraph::ArcIt it(digraph); it != INVALID; ++it) {
1528 1528
      if (digraph.source(it) == digraph.target(it)) return false;
1529 1529
    }
1530 1530
    return true;
1531 1531
  }
1532 1532

	
1533 1533
  /// \brief Returns true when there are not parallel edges in the graph.
1534 1534
  ///
1535 1535
  /// Returns true when there are not parallel edges in the graph.
1536 1536
  template <typename Digraph>
1537 1537
  bool parallelFree(const Digraph& digraph) {
1538 1538
    typename Digraph::template NodeMap<bool> reached(digraph, false);
1539 1539
    for (typename Digraph::NodeIt n(digraph); n != INVALID; ++n) {
1540 1540
      for (typename Digraph::OutArcIt a(digraph, n); a != INVALID; ++a) {
1541 1541
        if (reached[digraph.target(a)]) return false;
1542 1542
        reached.set(digraph.target(a), true);
1543 1543
      }
1544 1544
      for (typename Digraph::OutArcIt a(digraph, n); a != INVALID; ++a) {
1545 1545
        reached.set(digraph.target(a), false);
1546 1546
      }
1547 1547
    }
1548 1548
    return true;
1549 1549
  }
1550 1550

	
1551 1551
  /// \brief Returns true when there are not loop edges and parallel
1552 1552
  /// edges in the graph.
1553 1553
  ///
1554 1554
  /// Returns true when there are not loop edges and parallel edges in
1555 1555
  /// the graph.
1556 1556
  template <typename Digraph>
1557 1557
  bool simpleDigraph(const Digraph& digraph) {
1558 1558
    typename Digraph::template NodeMap<bool> reached(digraph, false);
1559 1559
    for (typename Digraph::NodeIt n(digraph); n != INVALID; ++n) {
1560 1560
      reached.set(n, true);
1561 1561
      for (typename Digraph::OutArcIt a(digraph, n); a != INVALID; ++a) {
1562 1562
        if (reached[digraph.target(a)]) return false;
1563 1563
        reached.set(digraph.target(a), true);
1564 1564
      }
1565 1565
      for (typename Digraph::OutArcIt a(digraph, n); a != INVALID; ++a) {
1566 1566
        reached.set(digraph.target(a), false);
1567 1567
      }
1568 1568
      reached.set(n, false);
1569 1569
    }
1570 1570
    return true;
1571 1571
  }
1572 1572

	
1573 1573
} //namespace lemon
1574 1574

	
1575 1575
#endif //LEMON_CONNECTIVITY_H
Ignore white space 1536 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2008
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_MAX_MATCHING_H
20 20
#define LEMON_MAX_MATCHING_H
21 21

	
22 22
#include <vector>
23 23
#include <queue>
24 24
#include <set>
25 25
#include <limits>
26 26

	
27 27
#include <lemon/core.h>
28 28
#include <lemon/unionfind.h>
29 29
#include <lemon/bin_heap.h>
30 30
#include <lemon/maps.h>
31 31

	
32 32
///\ingroup matching
33 33
///\file
34 34
///\brief Maximum matching algorithms in general graphs.
35 35

	
36 36
namespace lemon {
37 37

	
38 38
  /// \ingroup matching
39 39
  ///
40 40
  /// \brief Edmonds' alternating forest maximum matching algorithm.
41 41
  ///
42 42
  /// This class implements Edmonds' alternating forest matching
43 43
  /// algorithm. The algorithm can be started from an arbitrary initial
44 44
  /// matching (the default is the empty one)
45 45
  ///
46 46
  /// The dual solution of the problem is a map of the nodes to
47 47
  /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c
48 48
  /// MATCHED/C showing the Gallai-Edmonds decomposition of the
49 49
  /// graph. The nodes in \c EVEN/D induce a graph with
50 50
  /// factor-critical components, the nodes in \c ODD/A form the
51 51
  /// barrier, and the nodes in \c MATCHED/C induce a graph having a
52 52
  /// perfect matching. The number of the factor-critical components
53 53
  /// minus the number of barrier nodes is a lower bound on the
54 54
  /// unmatched nodes, and the matching is optimal if and only if this bound is
55 55
  /// tight. This decomposition can be attained by calling \c
56 56
  /// decomposition() after running the algorithm.
57 57
  ///
58 58
  /// \param _Graph The graph type the algorithm runs on.
59 59
  template <typename _Graph>
60 60
  class MaxMatching {
61 61
  public:
62 62

	
63 63
    typedef _Graph Graph;
64 64
    typedef typename Graph::template NodeMap<typename Graph::Arc>
65 65
    MatchingMap;
66 66

	
67 67
    ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
68 68
    ///
69 69
    ///Indicates the Gallai-Edmonds decomposition of the graph. The
70 70
    ///nodes with Status \c EVEN/D induce a graph with factor-critical
71 71
    ///components, the nodes in \c ODD/A form the canonical barrier,
72 72
    ///and the nodes in \c MATCHED/C induce a graph having a perfect
73 73
    ///matching.
74 74
    enum Status {
75 75
      EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2
76 76
    };
77 77

	
78 78
    typedef typename Graph::template NodeMap<Status> StatusMap;
79 79

	
80 80
  private:
81 81

	
82 82
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
83 83

	
84 84
    typedef UnionFindEnum<IntNodeMap> BlossomSet;
85 85
    typedef ExtendFindEnum<IntNodeMap> TreeSet;
86 86
    typedef RangeMap<Node> NodeIntMap;
87 87
    typedef MatchingMap EarMap;
88 88
    typedef std::vector<Node> NodeQueue;
89 89

	
90 90
    const Graph& _graph;
91 91
    MatchingMap* _matching;
92 92
    StatusMap* _status;
93 93

	
94 94
    EarMap* _ear;
95 95

	
96 96
    IntNodeMap* _blossom_set_index;
97 97
    BlossomSet* _blossom_set;
98 98
    NodeIntMap* _blossom_rep;
99 99

	
100 100
    IntNodeMap* _tree_set_index;
101 101
    TreeSet* _tree_set;
102 102

	
103 103
    NodeQueue _node_queue;
104 104
    int _process, _postpone, _last;
105 105

	
106 106
    int _node_num;
107 107

	
108 108
  private:
109 109

	
110 110
    void createStructures() {
111 111
      _node_num = countNodes(_graph);
112 112
      if (!_matching) {
113 113
        _matching = new MatchingMap(_graph);
114 114
      }
115 115
      if (!_status) {
116 116
        _status = new StatusMap(_graph);
117 117
      }
118 118
      if (!_ear) {
119 119
        _ear = new EarMap(_graph);
120 120
      }
121 121
      if (!_blossom_set) {
122 122
        _blossom_set_index = new IntNodeMap(_graph);
123 123
        _blossom_set = new BlossomSet(*_blossom_set_index);
124 124
      }
125 125
      if (!_blossom_rep) {
126 126
        _blossom_rep = new NodeIntMap(_node_num);
127 127
      }
128 128
      if (!_tree_set) {
129 129
        _tree_set_index = new IntNodeMap(_graph);
130 130
        _tree_set = new TreeSet(*_tree_set_index);
131 131
      }
132 132
      _node_queue.resize(_node_num);
133 133
    }
134 134

	
135 135
    void destroyStructures() {
136 136
      if (_matching) {
137 137
        delete _matching;
138 138
      }
139 139
      if (_status) {
140 140
        delete _status;
141 141
      }
142 142
      if (_ear) {
143 143
        delete _ear;
144 144
      }
145 145
      if (_blossom_set) {
146 146
        delete _blossom_set;
147 147
        delete _blossom_set_index;
148 148
      }
149 149
      if (_blossom_rep) {
150 150
        delete _blossom_rep;
151 151
      }
152 152
      if (_tree_set) {
153 153
        delete _tree_set_index;
154 154
        delete _tree_set;
155 155
      }
156 156
    }
157 157

	
158 158
    void processDense(const Node& n) {
159 159
      _process = _postpone = _last = 0;
160 160
      _node_queue[_last++] = n;
161 161

	
162 162
      while (_process != _last) {
163 163
        Node u = _node_queue[_process++];
164 164
        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
165 165
          Node v = _graph.target(a);
166 166
          if ((*_status)[v] == MATCHED) {
167 167
            extendOnArc(a);
168 168
          } else if ((*_status)[v] == UNMATCHED) {
169 169
            augmentOnArc(a);
170 170
            return;
171 171
          }
172 172
        }
173 173
      }
174 174

	
175 175
      while (_postpone != _last) {
176 176
        Node u = _node_queue[_postpone++];
177 177

	
178 178
        for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
179 179
          Node v = _graph.target(a);
180 180

	
181 181
          if ((*_status)[v] == EVEN) {
182 182
            if (_blossom_set->find(u) != _blossom_set->find(v)) {
183 183
              shrinkOnEdge(a);
184 184
            }
185 185
          }
186 186

	
187 187
          while (_process != _last) {
188 188
            Node w = _node_queue[_process++];
189 189
            for (OutArcIt b(_graph, w); b != INVALID; ++b) {
190 190
              Node x = _graph.target(b);
191 191
              if ((*_status)[x] == MATCHED) {
192 192
                extendOnArc(b);
193 193
              } else if ((*_status)[x] == UNMATCHED) {
194 194
                augmentOnArc(b);
195 195
                return;
196 196
              }
197 197
            }
198 198
          }
199 199
        }
200 200
      }
201 201
    }
202 202

	
203 203
    void processSparse(const Node& n) {
204 204
      _process = _last = 0;
205 205
      _node_queue[_last++] = n;
206 206
      while (_process != _last) {
207 207
        Node u = _node_queue[_process++];
208 208
        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
209 209
          Node v = _graph.target(a);
210 210

	
211 211
          if ((*_status)[v] == EVEN) {
212 212
            if (_blossom_set->find(u) != _blossom_set->find(v)) {
213 213
              shrinkOnEdge(a);
214 214
            }
215 215
          } else if ((*_status)[v] == MATCHED) {
216 216
            extendOnArc(a);
217 217
          } else if ((*_status)[v] == UNMATCHED) {
218 218
            augmentOnArc(a);
219 219
            return;
220 220
          }
221 221
        }
222 222
      }
223 223
    }
224 224

	
225 225
    void shrinkOnEdge(const Edge& e) {
226 226
      Node nca = INVALID;
227 227

	
228 228
      {
229 229
        std::set<Node> left_set, right_set;
230 230

	
231 231
        Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
232 232
        left_set.insert(left);
233 233

	
234 234
        Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
235 235
        right_set.insert(right);
236 236

	
237 237
        while (true) {
238 238
          if ((*_matching)[left] == INVALID) break;
239 239
          left = _graph.target((*_matching)[left]);
240 240
          left = (*_blossom_rep)[_blossom_set->
241 241
                                 find(_graph.target((*_ear)[left]))];
242 242
          if (right_set.find(left) != right_set.end()) {
243 243
            nca = left;
244 244
            break;
245 245
          }
246 246
          left_set.insert(left);
247 247

	
248 248
          if ((*_matching)[right] == INVALID) break;
249 249
          right = _graph.target((*_matching)[right]);
250 250
          right = (*_blossom_rep)[_blossom_set->
251 251
                                  find(_graph.target((*_ear)[right]))];
252 252
          if (left_set.find(right) != left_set.end()) {
253 253
            nca = right;
254 254
            break;
255 255
          }
256 256
          right_set.insert(right);
257 257
        }
258 258

	
259 259
        if (nca == INVALID) {
260 260
          if ((*_matching)[left] == INVALID) {
261 261
            nca = right;
262 262
            while (left_set.find(nca) == left_set.end()) {
263 263
              nca = _graph.target((*_matching)[nca]);
264 264
              nca =(*_blossom_rep)[_blossom_set->
265 265
                                   find(_graph.target((*_ear)[nca]))];
266 266
            }
267 267
          } else {
268 268
            nca = left;
269 269
            while (right_set.find(nca) == right_set.end()) {
270 270
              nca = _graph.target((*_matching)[nca]);
271 271
              nca = (*_blossom_rep)[_blossom_set->
272 272
                                   find(_graph.target((*_ear)[nca]))];
273 273
            }
274 274
          }
275 275
        }
276 276
      }
277 277

	
278 278
      {
279 279

	
280 280
        Node node = _graph.u(e);
281 281
        Arc arc = _graph.direct(e, true);
282 282
        Node base = (*_blossom_rep)[_blossom_set->find(node)];
283 283

	
284 284
        while (base != nca) {
285 285
          _ear->set(node, arc);
286 286

	
287 287
          Node n = node;
288 288
          while (n != base) {
289 289
            n = _graph.target((*_matching)[n]);
290 290
            Arc a = (*_ear)[n];
291 291
            n = _graph.target(a);
292 292
            _ear->set(n, _graph.oppositeArc(a));
293 293
          }
294 294
          node = _graph.target((*_matching)[base]);
295 295
          _tree_set->erase(base);
296 296
          _tree_set->erase(node);
297 297
          _blossom_set->insert(node, _blossom_set->find(base));
298 298
          _status->set(node, EVEN);
299 299
          _node_queue[_last++] = node;
300 300
          arc = _graph.oppositeArc((*_ear)[node]);
301 301
          node = _graph.target((*_ear)[node]);
302 302
          base = (*_blossom_rep)[_blossom_set->find(node)];
303 303
          _blossom_set->join(_graph.target(arc), base);
304 304
        }
305 305
      }
306 306

	
307 307
      _blossom_rep->set(_blossom_set->find(nca), nca);
308 308

	
309 309
      {
310 310

	
311 311
        Node node = _graph.v(e);
312 312
        Arc arc = _graph.direct(e, false);
313 313
        Node base = (*_blossom_rep)[_blossom_set->find(node)];
314 314

	
315 315
        while (base != nca) {
316 316
          _ear->set(node, arc);
317 317

	
318 318
          Node n = node;
319 319
          while (n != base) {
320 320
            n = _graph.target((*_matching)[n]);
321 321
            Arc a = (*_ear)[n];
322 322
            n = _graph.target(a);
323 323
            _ear->set(n, _graph.oppositeArc(a));
324 324
          }
325 325
          node = _graph.target((*_matching)[base]);
326 326
          _tree_set->erase(base);
327 327
          _tree_set->erase(node);
328 328
          _blossom_set->insert(node, _blossom_set->find(base));
329 329
          _status->set(node, EVEN);
330 330
          _node_queue[_last++] = node;
331 331
          arc = _graph.oppositeArc((*_ear)[node]);
332 332
          node = _graph.target((*_ear)[node]);
333 333
          base = (*_blossom_rep)[_blossom_set->find(node)];
334 334
          _blossom_set->join(_graph.target(arc), base);
335 335
        }
336 336
      }
337 337

	
338 338
      _blossom_rep->set(_blossom_set->find(nca), nca);
339 339
    }
340 340

	
341 341

	
342 342

	
343 343
    void extendOnArc(const Arc& a) {
344 344
      Node base = _graph.source(a);
345 345
      Node odd = _graph.target(a);
346 346

	
347 347
      _ear->set(odd, _graph.oppositeArc(a));
348 348
      Node even = _graph.target((*_matching)[odd]);
349 349
      _blossom_rep->set(_blossom_set->insert(even), even);
350 350
      _status->set(odd, ODD);
351 351
      _status->set(even, EVEN);
352 352
      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
353 353
      _tree_set->insert(odd, tree);
354 354
      _tree_set->insert(even, tree);
355 355
      _node_queue[_last++] = even;
356 356

	
357 357
    }
358 358

	
359 359
    void augmentOnArc(const Arc& a) {
360 360
      Node even = _graph.source(a);
361 361
      Node odd = _graph.target(a);
362 362

	
363 363
      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
364 364

	
365 365
      _matching->set(odd, _graph.oppositeArc(a));
366 366
      _status->set(odd, MATCHED);
367 367

	
368 368
      Arc arc = (*_matching)[even];
369 369
      _matching->set(even, a);
370 370

	
371 371
      while (arc != INVALID) {
372 372
        odd = _graph.target(arc);
373 373
        arc = (*_ear)[odd];
374 374
        even = _graph.target(arc);
375 375
        _matching->set(odd, arc);
376 376
        arc = (*_matching)[even];
377 377
        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
378 378
      }
379 379

	
380 380
      for (typename TreeSet::ItemIt it(*_tree_set, tree);
381 381
           it != INVALID; ++it) {
382 382
        if ((*_status)[it] == ODD) {
383 383
          _status->set(it, MATCHED);
384 384
        } else {
385 385
          int blossom = _blossom_set->find(it);
386 386
          for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
387 387
               jt != INVALID; ++jt) {
388 388
            _status->set(jt, MATCHED);
389 389
          }
390 390
          _blossom_set->eraseClass(blossom);
391 391
        }
392 392
      }
393 393
      _tree_set->eraseClass(tree);
394 394

	
395 395
    }
396 396

	
397 397
  public:
398 398

	
399 399
    /// \brief Constructor
400 400
    ///
401 401
    /// Constructor.
402 402
    MaxMatching(const Graph& graph)
403 403
      : _graph(graph), _matching(0), _status(0), _ear(0),
404 404
        _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
405 405
        _tree_set_index(0), _tree_set(0) {}
406 406

	
407 407
    ~MaxMatching() {
408 408
      destroyStructures();
409 409
    }
410 410

	
411 411
    /// \name Execution control
412 412
    /// The simplest way to execute the algorithm is to use the
413 413
    /// \c run() member function.
414 414
    /// \n
415 415

	
416 416
    /// If you need better control on the execution, you must call
417 417
    /// \ref init(), \ref greedyInit() or \ref matchingInit()
418 418
    /// functions first, then you can start the algorithm with the \ref
419
    /// startParse() or startDense() functions.
419
    /// startSparse() or startDense() functions.
420 420

	
421 421
    ///@{
422 422

	
423 423
    /// \brief Sets the actual matching to the empty matching.
424 424
    ///
425 425
    /// Sets the actual matching to the empty matching.
426 426
    ///
427 427
    void init() {
428 428
      createStructures();
429 429
      for(NodeIt n(_graph); n != INVALID; ++n) {
430 430
        _matching->set(n, INVALID);
431 431
        _status->set(n, UNMATCHED);
432 432
      }
433 433
    }
434 434

	
435 435
    ///\brief Finds an initial matching in a greedy way
436 436
    ///
437 437
    ///It finds an initial matching in a greedy way.
438 438
    void greedyInit() {
439 439
      createStructures();
440 440
      for (NodeIt n(_graph); n != INVALID; ++n) {
441 441
        _matching->set(n, INVALID);
442 442
        _status->set(n, UNMATCHED);
443 443
      }
444 444
      for (NodeIt n(_graph); n != INVALID; ++n) {
445 445
        if ((*_matching)[n] == INVALID) {
446 446
          for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
447 447
            Node v = _graph.target(a);
448 448
            if ((*_matching)[v] == INVALID && v != n) {
449 449
              _matching->set(n, a);
450 450
              _status->set(n, MATCHED);
451 451
              _matching->set(v, _graph.oppositeArc(a));
452 452
              _status->set(v, MATCHED);
453 453
              break;
454 454
            }
455 455
          }
456 456
        }
457 457
      }
458 458
    }
459 459

	
460 460

	
461 461
    /// \brief Initialize the matching from a map containing.
462 462
    ///
463 463
    /// Initialize the matching from a \c bool valued \c Edge map. This
464 464
    /// map must have the property that there are no two incident edges
465 465
    /// with true value, ie. it contains a matching.
466 466
    /// \return %True if the map contains a matching.
467 467
    template <typename MatchingMap>
468 468
    bool matchingInit(const MatchingMap& matching) {
469 469
      createStructures();
470 470

	
471 471
      for (NodeIt n(_graph); n != INVALID; ++n) {
472 472
        _matching->set(n, INVALID);
473 473
        _status->set(n, UNMATCHED);
474 474
      }
475 475
      for(EdgeIt e(_graph); e!=INVALID; ++e) {
476 476
        if (matching[e]) {
477 477

	
478 478
          Node u = _graph.u(e);
479 479
          if ((*_matching)[u] != INVALID) return false;
480 480
          _matching->set(u, _graph.direct(e, true));
481 481
          _status->set(u, MATCHED);
482 482

	
483 483
          Node v = _graph.v(e);
484 484
          if ((*_matching)[v] != INVALID) return false;
485 485
          _matching->set(v, _graph.direct(e, false));
486 486
          _status->set(v, MATCHED);
487 487
        }
488 488
      }
489 489
      return true;
490 490
    }
491 491

	
492 492
    /// \brief Starts Edmonds' algorithm
493 493
    ///
494 494
    /// If runs the original Edmonds' algorithm.
495 495
    void startSparse() {
496 496
      for(NodeIt n(_graph); n != INVALID; ++n) {
497 497
        if ((*_status)[n] == UNMATCHED) {
498 498
          (*_blossom_rep)[_blossom_set->insert(n)] = n;
499 499
          _tree_set->insert(n);
500 500
          _status->set(n, EVEN);
501 501
          processSparse(n);
502 502
        }
503 503
      }
504 504
    }
505 505

	
506 506
    /// \brief Starts Edmonds' algorithm.
507 507
    ///
508 508
    /// It runs Edmonds' algorithm with a heuristic of postponing
509 509
    /// shrinks, therefore resulting in a faster algorithm for dense graphs.
510 510
    void startDense() {
511 511
      for(NodeIt n(_graph); n != INVALID; ++n) {
512 512
        if ((*_status)[n] == UNMATCHED) {
513 513
          (*_blossom_rep)[_blossom_set->insert(n)] = n;
514 514
          _tree_set->insert(n);
515 515
          _status->set(n, EVEN);
516 516
          processDense(n);
517 517
        }
518 518
      }
519 519
    }
520 520

	
521 521

	
522 522
    /// \brief Runs Edmonds' algorithm
523 523
    ///
524 524
    /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>)
525 525
    /// or Edmonds' algorithm with a heuristic of
526 526
    /// postponing shrinks for dense graphs.
527 527
    void run() {
528 528
      if (countEdges(_graph) < 2 * countNodes(_graph)) {
529 529
        greedyInit();
530 530
        startSparse();
531 531
      } else {
532 532
        init();
533 533
        startDense();
534 534
      }
535 535
    }
536 536

	
537 537
    /// @}
538 538

	
539 539
    /// \name Primal solution
540 540
    /// Functions to get the primal solution, ie. the matching.
541 541

	
542 542
    /// @{
543 543

	
544 544
    ///\brief Returns the size of the current matching.
545 545
    ///
546 546
    ///Returns the size of the current matching. After \ref
547 547
    ///run() it returns the size of the maximum matching in the graph.
548 548
    int matchingSize() const {
549 549
      int size = 0;
550 550
      for (NodeIt n(_graph); n != INVALID; ++n) {
551 551
        if ((*_matching)[n] != INVALID) {
552 552
          ++size;
553 553
        }
554 554
      }
555 555
      return size / 2;
556 556
    }
557 557

	
558 558
    /// \brief Returns true when the edge is in the matching.
559 559
    ///
560 560
    /// Returns true when the edge is in the matching.
561 561
    bool matching(const Edge& edge) const {
562 562
      return edge == (*_matching)[_graph.u(edge)];
563 563
    }
564 564

	
565 565
    /// \brief Returns the matching edge incident to the given node.
566 566
    ///
567 567
    /// Returns the matching edge of a \c node in the actual matching or
568 568
    /// INVALID if the \c node is not covered by the actual matching.
569 569
    Arc matching(const Node& n) const {
570 570
      return (*_matching)[n];
571 571
    }
572 572

	
573 573
    ///\brief Returns the mate of a node in the actual matching.
574 574
    ///
575 575
    ///Returns the mate of a \c node in the actual matching or
576 576
    ///INVALID if the \c node is not covered by the actual matching.
577 577
    Node mate(const Node& n) const {
578 578
      return (*_matching)[n] != INVALID ?
579 579
        _graph.target((*_matching)[n]) : INVALID;
580 580
    }
581 581

	
582 582
    /// @}
583 583

	
584 584
    /// \name Dual solution
585 585
    /// Functions to get the dual solution, ie. the decomposition.
586 586

	
587 587
    /// @{
588 588

	
589 589
    /// \brief Returns the class of the node in the Edmonds-Gallai
590 590
    /// decomposition.
591 591
    ///
592 592
    /// Returns the class of the node in the Edmonds-Gallai
593 593
    /// decomposition.
594 594
    Status decomposition(const Node& n) const {
595 595
      return (*_status)[n];
596 596
    }
597 597

	
598 598
    /// \brief Returns true when the node is in the barrier.
599 599
    ///
600 600
    /// Returns true when the node is in the barrier.
601 601
    bool barrier(const Node& n) const {
602 602
      return (*_status)[n] == ODD;
603 603
    }
604 604

	
605 605
    /// @}
606 606

	
607 607
  };
608 608

	
609 609
  /// \ingroup matching
610 610
  ///
611 611
  /// \brief Weighted matching in general graphs
612 612
  ///
613 613
  /// This class provides an efficient implementation of Edmond's
614 614
  /// maximum weighted matching algorithm. The implementation is based
615 615
  /// on extensive use of priority queues and provides
616 616
  /// \f$O(nm\log(n))\f$ time complexity.
617 617
  ///
618 618
  /// The maximum weighted matching problem is to find undirected
619 619
  /// edges in the graph with maximum overall weight and no two of
620 620
  /// them shares their ends. The problem can be formulated with the
621 621
  /// following linear program.
622 622
  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
623 623
  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
624 624
      \quad \forall B\in\mathcal{O}\f] */
625 625
  /// \f[x_e \ge 0\quad \forall e\in E\f]
626 626
  /// \f[\max \sum_{e\in E}x_ew_e\f]
627 627
  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
628 628
  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
629 629
  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
630 630
  /// subsets of the nodes.
631 631
  ///
632 632
  /// The algorithm calculates an optimal matching and a proof of the
633 633
  /// optimality. The solution of the dual problem can be used to check
634 634
  /// the result of the algorithm. The dual linear problem is the
635 635
  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
636 636
      z_B \ge w_{uv} \quad \forall uv\in E\f] */
637 637
  /// \f[y_u \ge 0 \quad \forall u \in V\f]
638 638
  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
639 639
  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
640 640
      \frac{\vert B \vert - 1}{2}z_B\f] */
641 641
  ///
642 642
  /// The algorithm can be executed with \c run() or the \c init() and
643 643
  /// then the \c start() member functions. After it the matching can
644 644
  /// be asked with \c matching() or mate() functions. The dual
645 645
  /// solution can be get with \c nodeValue(), \c blossomNum() and \c
646 646
  /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
647 647
  /// "BlossomIt" nested class, which is able to iterate on the nodes
648 648
  /// of a blossom. If the value type is integral then the dual
649 649
  /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
650 650
  template <typename _Graph,
651 651
            typename _WeightMap = typename _Graph::template EdgeMap<int> >
652 652
  class MaxWeightedMatching {
653 653
  public:
654 654

	
655 655
    typedef _Graph Graph;
656 656
    typedef _WeightMap WeightMap;
657 657
    typedef typename WeightMap::Value Value;
658 658

	
659 659
    /// \brief Scaling factor for dual solution
660 660
    ///
661 661
    /// Scaling factor for dual solution, it is equal to 4 or 1
662 662
    /// according to the value type.
663 663
    static const int dualScale =
664 664
      std::numeric_limits<Value>::is_integer ? 4 : 1;
665 665

	
666 666
    typedef typename Graph::template NodeMap<typename Graph::Arc>
667 667
    MatchingMap;
668 668

	
669 669
  private:
670 670

	
671 671
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
672 672

	
673 673
    typedef typename Graph::template NodeMap<Value> NodePotential;
674 674
    typedef std::vector<Node> BlossomNodeList;
675 675

	
676 676
    struct BlossomVariable {
677 677
      int begin, end;
678 678
      Value value;
679 679

	
680 680
      BlossomVariable(int _begin, int _end, Value _value)
681 681
        : begin(_begin), end(_end), value(_value) {}
682 682

	
683 683
    };
684 684

	
685 685
    typedef std::vector<BlossomVariable> BlossomPotential;
686 686

	
687 687
    const Graph& _graph;
688 688
    const WeightMap& _weight;
689 689

	
690 690
    MatchingMap* _matching;
691 691

	
692 692
    NodePotential* _node_potential;
693 693

	
694 694
    BlossomPotential _blossom_potential;
695 695
    BlossomNodeList _blossom_node_list;
696 696

	
697 697
    int _node_num;
698 698
    int _blossom_num;
699 699

	
700 700
    typedef RangeMap<int> IntIntMap;
701 701

	
702 702
    enum Status {
703 703
      EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
704 704
    };
705 705

	
706 706
    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
707 707
    struct BlossomData {
708 708
      int tree;
709 709
      Status status;
710 710
      Arc pred, next;
711 711
      Value pot, offset;
712 712
      Node base;
713 713
    };
714 714

	
715 715
    IntNodeMap *_blossom_index;
716 716
    BlossomSet *_blossom_set;
717 717
    RangeMap<BlossomData>* _blossom_data;
718 718

	
719 719
    IntNodeMap *_node_index;
720 720
    IntArcMap *_node_heap_index;
721 721

	
722 722
    struct NodeData {
723 723

	
724 724
      NodeData(IntArcMap& node_heap_index)
725 725
        : heap(node_heap_index) {}
726 726

	
727 727
      int blossom;
728 728
      Value pot;
729 729
      BinHeap<Value, IntArcMap> heap;
730 730
      std::map<int, Arc> heap_index;
731 731

	
732 732
      int tree;
733 733
    };
734 734

	
735 735
    RangeMap<NodeData>* _node_data;
736 736

	
737 737
    typedef ExtendFindEnum<IntIntMap> TreeSet;
738 738

	
739 739
    IntIntMap *_tree_set_index;
740 740
    TreeSet *_tree_set;
741 741

	
742 742
    IntNodeMap *_delta1_index;
743 743
    BinHeap<Value, IntNodeMap> *_delta1;
744 744

	
745 745
    IntIntMap *_delta2_index;
746 746
    BinHeap<Value, IntIntMap> *_delta2;
747 747

	
748 748
    IntEdgeMap *_delta3_index;
749 749
    BinHeap<Value, IntEdgeMap> *_delta3;
750 750

	
751 751
    IntIntMap *_delta4_index;
752 752
    BinHeap<Value, IntIntMap> *_delta4;
753 753

	
754 754
    Value _delta_sum;
755 755

	
756 756
    void createStructures() {
757 757
      _node_num = countNodes(_graph);
758 758
      _blossom_num = _node_num * 3 / 2;
759 759

	
760 760
      if (!_matching) {
761 761
        _matching = new MatchingMap(_graph);
762 762
      }
763 763
      if (!_node_potential) {
764 764
        _node_potential = new NodePotential(_graph);
765 765
      }
766 766
      if (!_blossom_set) {
767 767
        _blossom_index = new IntNodeMap(_graph);
768 768
        _blossom_set = new BlossomSet(*_blossom_index);
769 769
        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
770 770
      }
771 771

	
772 772
      if (!_node_index) {
773 773
        _node_index = new IntNodeMap(_graph);
774 774
        _node_heap_index = new IntArcMap(_graph);
775 775
        _node_data = new RangeMap<NodeData>(_node_num,
776 776
                                              NodeData(*_node_heap_index));
777 777
      }
778 778

	
779 779
      if (!_tree_set) {
780 780
        _tree_set_index = new IntIntMap(_blossom_num);
781 781
        _tree_set = new TreeSet(*_tree_set_index);
782 782
      }
783 783
      if (!_delta1) {
784 784
        _delta1_index = new IntNodeMap(_graph);
785 785
        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
786 786
      }
787 787
      if (!_delta2) {
788 788
        _delta2_index = new IntIntMap(_blossom_num);
789 789
        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
790 790
      }
791 791
      if (!_delta3) {
792 792
        _delta3_index = new IntEdgeMap(_graph);
793 793
        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
794 794
      }
795 795
      if (!_delta4) {
796 796
        _delta4_index = new IntIntMap(_blossom_num);
797 797
        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
798 798
      }
799 799
    }
800 800

	
801 801
    void destroyStructures() {
802 802
      _node_num = countNodes(_graph);
803 803
      _blossom_num = _node_num * 3 / 2;
804 804

	
805 805
      if (_matching) {
806 806
        delete _matching;
807 807
      }
808 808
      if (_node_potential) {
809 809
        delete _node_potential;
810 810
      }
811 811
      if (_blossom_set) {
812 812
        delete _blossom_index;
813 813
        delete _blossom_set;
814 814
        delete _blossom_data;
815 815
      }
816 816

	
817 817
      if (_node_index) {
818 818
        delete _node_index;
819 819
        delete _node_heap_index;
820 820
        delete _node_data;
821 821
      }
822 822

	
823 823
      if (_tree_set) {
824 824
        delete _tree_set_index;
825 825
        delete _tree_set;
826 826
      }
827 827
      if (_delta1) {
828 828
        delete _delta1_index;
829 829
        delete _delta1;
830 830
      }
831 831
      if (_delta2) {
832 832
        delete _delta2_index;
833 833
        delete _delta2;
834 834
      }
835 835
      if (_delta3) {
836 836
        delete _delta3_index;
837 837
        delete _delta3;
838 838
      }
839 839
      if (_delta4) {
840 840
        delete _delta4_index;
841 841
        delete _delta4;
842 842
      }
843 843
    }
844 844

	
845 845
    void matchedToEven(int blossom, int tree) {
846 846
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
847 847
        _delta2->erase(blossom);
848 848
      }
849 849

	
850 850
      if (!_blossom_set->trivial(blossom)) {
851 851
        (*_blossom_data)[blossom].pot -=
852 852
          2 * (_delta_sum - (*_blossom_data)[blossom].offset);
853 853
      }
854 854

	
855 855
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
856 856
           n != INVALID; ++n) {
857 857

	
858 858
        _blossom_set->increase(n, std::numeric_limits<Value>::max());
859 859
        int ni = (*_node_index)[n];
860 860

	
861 861
        (*_node_data)[ni].heap.clear();
862 862
        (*_node_data)[ni].heap_index.clear();
863 863

	
864 864
        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
865 865

	
866 866
        _delta1->push(n, (*_node_data)[ni].pot);
867 867

	
868 868
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
869 869
          Node v = _graph.source(e);
870 870
          int vb = _blossom_set->find(v);
871 871
          int vi = (*_node_index)[v];
872 872

	
873 873
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
874 874
            dualScale * _weight[e];
875 875

	
876 876
          if ((*_blossom_data)[vb].status == EVEN) {
877 877
            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
878 878
              _delta3->push(e, rw / 2);
879 879
            }
880 880
          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
881 881
            if (_delta3->state(e) != _delta3->IN_HEAP) {
882 882
              _delta3->push(e, rw);
883 883
            }
884 884
          } else {
885 885
            typename std::map<int, Arc>::iterator it =
886 886
              (*_node_data)[vi].heap_index.find(tree);
887 887

	
888 888
            if (it != (*_node_data)[vi].heap_index.end()) {
889 889
              if ((*_node_data)[vi].heap[it->second] > rw) {
890 890
                (*_node_data)[vi].heap.replace(it->second, e);
891 891
                (*_node_data)[vi].heap.decrease(e, rw);
892 892
                it->second = e;
893 893
              }
894 894
            } else {
895 895
              (*_node_data)[vi].heap.push(e, rw);
896 896
              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
897 897
            }
898 898

	
899 899
            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
900 900
              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
901 901

	
902 902
              if ((*_blossom_data)[vb].status == MATCHED) {
903 903
                if (_delta2->state(vb) != _delta2->IN_HEAP) {
904 904
                  _delta2->push(vb, _blossom_set->classPrio(vb) -
905 905
                               (*_blossom_data)[vb].offset);
906 906
                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
907 907
                           (*_blossom_data)[vb].offset){
908 908
                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
909 909
                                   (*_blossom_data)[vb].offset);
910 910
                }
911 911
              }
912 912
            }
913 913
          }
914 914
        }
915 915
      }
916 916
      (*_blossom_data)[blossom].offset = 0;
917 917
    }
918 918

	
919 919
    void matchedToOdd(int blossom) {
920 920
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
921 921
        _delta2->erase(blossom);
922 922
      }
923 923
      (*_blossom_data)[blossom].offset += _delta_sum;
924 924
      if (!_blossom_set->trivial(blossom)) {
925 925
        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
926 926
                     (*_blossom_data)[blossom].offset);
927 927
      }
928 928
    }
929 929

	
930 930
    void evenToMatched(int blossom, int tree) {
931 931
      if (!_blossom_set->trivial(blossom)) {
932 932
        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
933 933
      }
934 934

	
935 935
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
936 936
           n != INVALID; ++n) {
937 937
        int ni = (*_node_index)[n];
938 938
        (*_node_data)[ni].pot -= _delta_sum;
939 939

	
940 940
        _delta1->erase(n);
941 941

	
942 942
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
943 943
          Node v = _graph.source(e);
944 944
          int vb = _blossom_set->find(v);
945 945
          int vi = (*_node_index)[v];
946 946

	
947 947
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
948 948
            dualScale * _weight[e];
949 949

	
950 950
          if (vb == blossom) {
951 951
            if (_delta3->state(e) == _delta3->IN_HEAP) {
952 952
              _delta3->erase(e);
953 953
            }
954 954
          } else if ((*_blossom_data)[vb].status == EVEN) {
955 955

	
956 956
            if (_delta3->state(e) == _delta3->IN_HEAP) {
957 957
              _delta3->erase(e);
958 958
            }
959 959

	
960 960
            int vt = _tree_set->find(vb);
961 961

	
962 962
            if (vt != tree) {
963 963

	
964 964
              Arc r = _graph.oppositeArc(e);
965 965

	
966 966
              typename std::map<int, Arc>::iterator it =
967 967
                (*_node_data)[ni].heap_index.find(vt);
968 968

	
969 969
              if (it != (*_node_data)[ni].heap_index.end()) {
970 970
                if ((*_node_data)[ni].heap[it->second] > rw) {
971 971
                  (*_node_data)[ni].heap.replace(it->second, r);
972 972
                  (*_node_data)[ni].heap.decrease(r, rw);
973 973
                  it->second = r;
974 974
                }
975 975
              } else {
976 976
                (*_node_data)[ni].heap.push(r, rw);
977 977
                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
978 978
              }
979 979

	
980 980
              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
981 981
                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
982 982

	
983 983
                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
984 984
                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
985 985
                               (*_blossom_data)[blossom].offset);
986 986
                } else if ((*_delta2)[blossom] >
987 987
                           _blossom_set->classPrio(blossom) -
988 988
                           (*_blossom_data)[blossom].offset){
989 989
                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
990 990
                                   (*_blossom_data)[blossom].offset);
991 991
                }
992 992
              }
993 993
            }
994 994

	
995 995
          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
996 996
            if (_delta3->state(e) == _delta3->IN_HEAP) {
997 997
              _delta3->erase(e);
998 998
            }
999 999
          } else {
1000 1000

	
1001 1001
            typename std::map<int, Arc>::iterator it =
1002 1002
              (*_node_data)[vi].heap_index.find(tree);
1003 1003

	
1004 1004
            if (it != (*_node_data)[vi].heap_index.end()) {
1005 1005
              (*_node_data)[vi].heap.erase(it->second);
1006 1006
              (*_node_data)[vi].heap_index.erase(it);
1007 1007
              if ((*_node_data)[vi].heap.empty()) {
1008 1008
                _blossom_set->increase(v, std::numeric_limits<Value>::max());
1009 1009
              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1010 1010
                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1011 1011
              }
1012 1012

	
1013 1013
              if ((*_blossom_data)[vb].status == MATCHED) {
1014 1014
                if (_blossom_set->classPrio(vb) ==
1015 1015
                    std::numeric_limits<Value>::max()) {
1016 1016
                  _delta2->erase(vb);
1017 1017
                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1018 1018
                           (*_blossom_data)[vb].offset) {
1019 1019
                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
1020 1020
                                   (*_blossom_data)[vb].offset);
1021 1021
                }
1022 1022
              }
1023 1023
            }
1024 1024
          }
1025 1025
        }
1026 1026
      }
1027 1027
    }
1028 1028

	
1029 1029
    void oddToMatched(int blossom) {
1030 1030
      (*_blossom_data)[blossom].offset -= _delta_sum;
1031 1031

	
1032 1032
      if (_blossom_set->classPrio(blossom) !=
1033 1033
          std::numeric_limits<Value>::max()) {
1034 1034
        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1035 1035
                       (*_blossom_data)[blossom].offset);
1036 1036
      }
1037 1037

	
1038 1038
      if (!_blossom_set->trivial(blossom)) {
1039 1039
        _delta4->erase(blossom);
1040 1040
      }
1041 1041
    }
1042 1042

	
1043 1043
    void oddToEven(int blossom, int tree) {
1044 1044
      if (!_blossom_set->trivial(blossom)) {
1045 1045
        _delta4->erase(blossom);
1046 1046
        (*_blossom_data)[blossom].pot -=
1047 1047
          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1048 1048
      }
1049 1049

	
1050 1050
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1051 1051
           n != INVALID; ++n) {
1052 1052
        int ni = (*_node_index)[n];
1053 1053

	
1054 1054
        _blossom_set->increase(n, std::numeric_limits<Value>::max());
1055 1055

	
1056 1056
        (*_node_data)[ni].heap.clear();
1057 1057
        (*_node_data)[ni].heap_index.clear();
1058 1058
        (*_node_data)[ni].pot +=
1059 1059
          2 * _delta_sum - (*_blossom_data)[blossom].offset;
1060 1060

	
1061 1061
        _delta1->push(n, (*_node_data)[ni].pot);
1062 1062

	
1063 1063
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
1064 1064
          Node v = _graph.source(e);
1065 1065
          int vb = _blossom_set->find(v);
1066 1066
          int vi = (*_node_index)[v];
1067 1067

	
1068 1068
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1069 1069
            dualScale * _weight[e];
1070 1070

	
1071 1071
          if ((*_blossom_data)[vb].status == EVEN) {
1072 1072
            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1073 1073
              _delta3->push(e, rw / 2);
1074 1074
            }
1075 1075
          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1076 1076
            if (_delta3->state(e) != _delta3->IN_HEAP) {
1077 1077
              _delta3->push(e, rw);
1078 1078
            }
1079 1079
          } else {
1080 1080

	
1081 1081
            typename std::map<int, Arc>::iterator it =
1082 1082
              (*_node_data)[vi].heap_index.find(tree);
1083 1083

	
1084 1084
            if (it != (*_node_data)[vi].heap_index.end()) {
1085 1085
              if ((*_node_data)[vi].heap[it->second] > rw) {
1086 1086
                (*_node_data)[vi].heap.replace(it->second, e);
1087 1087
                (*_node_data)[vi].heap.decrease(e, rw);
1088 1088
                it->second = e;
1089 1089
              }
1090 1090
            } else {
1091 1091
              (*_node_data)[vi].heap.push(e, rw);
1092 1092
              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1093 1093
            }
1094 1094

	
1095 1095
            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1096 1096
              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1097 1097

	
1098 1098
              if ((*_blossom_data)[vb].status == MATCHED) {
1099 1099
                if (_delta2->state(vb) != _delta2->IN_HEAP) {
1100 1100
                  _delta2->push(vb, _blossom_set->classPrio(vb) -
1101 1101
                               (*_blossom_data)[vb].offset);
1102 1102
                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1103 1103
                           (*_blossom_data)[vb].offset) {
1104 1104
                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1105 1105
                                   (*_blossom_data)[vb].offset);
1106 1106
                }
1107 1107
              }
1108 1108
            }
1109 1109
          }
1110 1110
        }
1111 1111
      }
1112 1112
      (*_blossom_data)[blossom].offset = 0;
1113 1113
    }
1114 1114

	
1115 1115

	
1116 1116
    void matchedToUnmatched(int blossom) {
1117 1117
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1118 1118
        _delta2->erase(blossom);
1119 1119
      }
1120 1120

	
1121 1121
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1122 1122
           n != INVALID; ++n) {
1123 1123
        int ni = (*_node_index)[n];
1124 1124

	
1125 1125
        _blossom_set->increase(n, std::numeric_limits<Value>::max());
1126 1126

	
1127 1127
        (*_node_data)[ni].heap.clear();
1128 1128
        (*_node_data)[ni].heap_index.clear();
1129 1129

	
1130 1130
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1131 1131
          Node v = _graph.target(e);
1132 1132
          int vb = _blossom_set->find(v);
1133 1133
          int vi = (*_node_index)[v];
1134 1134

	
1135 1135
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1136 1136
            dualScale * _weight[e];
1137 1137

	
1138 1138
          if ((*_blossom_data)[vb].status == EVEN) {
1139 1139
            if (_delta3->state(e) != _delta3->IN_HEAP) {
1140 1140
              _delta3->push(e, rw);
1141 1141
            }
1142 1142
          }
1143 1143
        }
1144 1144
      }
1145 1145
    }
1146 1146

	
1147 1147
    void unmatchedToMatched(int blossom) {
1148 1148
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1149 1149
           n != INVALID; ++n) {
1150 1150
        int ni = (*_node_index)[n];
1151 1151

	
1152 1152
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
1153 1153
          Node v = _graph.source(e);
1154 1154
          int vb = _blossom_set->find(v);
1155 1155
          int vi = (*_node_index)[v];
1156 1156

	
1157 1157
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1158 1158
            dualScale * _weight[e];
1159 1159

	
1160 1160
          if (vb == blossom) {
1161 1161
            if (_delta3->state(e) == _delta3->IN_HEAP) {
1162 1162
              _delta3->erase(e);
1163 1163
            }
1164 1164
          } else if ((*_blossom_data)[vb].status == EVEN) {
1165 1165

	
1166 1166
            if (_delta3->state(e) == _delta3->IN_HEAP) {
1167 1167
              _delta3->erase(e);
1168 1168
            }
1169 1169

	
1170 1170
            int vt = _tree_set->find(vb);
1171 1171

	
1172 1172
            Arc r = _graph.oppositeArc(e);
1173 1173

	
1174 1174
            typename std::map<int, Arc>::iterator it =
1175 1175
              (*_node_data)[ni].heap_index.find(vt);
1176 1176

	
1177 1177
            if (it != (*_node_data)[ni].heap_index.end()) {
1178 1178
              if ((*_node_data)[ni].heap[it->second] > rw) {
1179 1179
                (*_node_data)[ni].heap.replace(it->second, r);
1180 1180
                (*_node_data)[ni].heap.decrease(r, rw);
1181 1181
                it->second = r;
1182 1182
              }
1183 1183
            } else {
1184 1184
              (*_node_data)[ni].heap.push(r, rw);
1185 1185
              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1186 1186
            }
1187 1187

	
Ignore white space 1536 line context
1 1
/* -*- C++ -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library
4 4
 *
5 5
 * Copyright (C) 2003-2008
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_SUURBALLE_H
20 20
#define LEMON_SUURBALLE_H
21 21

	
22 22
///\ingroup shortest_path
23 23
///\file
24 24
///\brief An algorithm for finding arc-disjoint paths between two
25 25
/// nodes having minimum total length.
26 26

	
27 27
#include <vector>
28 28
#include <lemon/bin_heap.h>
29 29
#include <lemon/path.h>
30 30

	
31 31
namespace lemon {
32 32

	
33 33
  /// \addtogroup shortest_path
34 34
  /// @{
35 35

	
36 36
  /// \brief Algorithm for finding arc-disjoint paths between two nodes
37 37
  /// having minimum total length.
38 38
  ///
39 39
  /// \ref lemon::Suurballe "Suurballe" implements an algorithm for
40 40
  /// finding arc-disjoint paths having minimum total length (cost)
41 41
  /// from a given source node to a given target node in a digraph.
42 42
  ///
43 43
  /// In fact, this implementation is the specialization of the
44 44
  /// \ref CapacityScaling "successive shortest path" algorithm.
45 45
  ///
46 46
  /// \tparam Digraph The digraph type the algorithm runs on.
47 47
  /// The default value is \c ListDigraph.
48 48
  /// \tparam LengthMap The type of the length (cost) map.
49 49
  /// The default value is <tt>Digraph::ArcMap<int></tt>.
50 50
  ///
51 51
  /// \warning Length values should be \e non-negative \e integers.
52 52
  ///
53 53
  /// \note For finding node-disjoint paths this algorithm can be used
54
  /// with \ref SplitDigraphAdaptor.
54
  /// with \ref SplitNodes.
55 55
#ifdef DOXYGEN
56 56
  template <typename Digraph, typename LengthMap>
57 57
#else
58 58
  template < typename Digraph = ListDigraph,
59 59
             typename LengthMap = typename Digraph::template ArcMap<int> >
60 60
#endif
61 61
  class Suurballe
62 62
  {
63 63
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
64 64

	
65 65
    typedef typename LengthMap::Value Length;
66 66
    typedef ConstMap<Arc, int> ConstArcMap;
67 67
    typedef typename Digraph::template NodeMap<Arc> PredMap;
68 68

	
69 69
  public:
70 70

	
71 71
    /// The type of the flow map.
72 72
    typedef typename Digraph::template ArcMap<int> FlowMap;
73 73
    /// The type of the potential map.
74 74
    typedef typename Digraph::template NodeMap<Length> PotentialMap;
75 75
    /// The type of the path structures.
76 76
    typedef SimplePath<Digraph> Path;
77 77

	
78 78
  private:
79 79
  
80 80
    /// \brief Special implementation of the Dijkstra algorithm
81 81
    /// for finding shortest paths in the residual network.
82 82
    ///
83 83
    /// \ref ResidualDijkstra is a special implementation of the
84 84
    /// \ref Dijkstra algorithm for finding shortest paths in the
85 85
    /// residual network of the digraph with respect to the reduced arc
86 86
    /// lengths and modifying the node potentials according to the
87 87
    /// distance of the nodes.
88 88
    class ResidualDijkstra
89 89
    {
90 90
      typedef typename Digraph::template NodeMap<int> HeapCrossRef;
91 91
      typedef BinHeap<Length, HeapCrossRef> Heap;
92 92

	
93 93
    private:
94 94

	
95 95
      // The digraph the algorithm runs on
96 96
      const Digraph &_graph;
97 97

	
98 98
      // The main maps
99 99
      const FlowMap &_flow;
100 100
      const LengthMap &_length;
101 101
      PotentialMap &_potential;
102 102

	
103 103
      // The distance map
104 104
      PotentialMap _dist;
105 105
      // The pred arc map
106 106
      PredMap &_pred;
107 107
      // The processed (i.e. permanently labeled) nodes
108 108
      std::vector<Node> _proc_nodes;
109 109
      
110 110
      Node _s;
111 111
      Node _t;
112 112

	
113 113
    public:
114 114

	
115 115
      /// Constructor.
116 116
      ResidualDijkstra( const Digraph &digraph,
117 117
                        const FlowMap &flow,
118 118
                        const LengthMap &length,
119 119
                        PotentialMap &potential,
120 120
                        PredMap &pred,
121 121
                        Node s, Node t ) :
122 122
        _graph(digraph), _flow(flow), _length(length), _potential(potential),
123 123
        _dist(digraph), _pred(pred), _s(s), _t(t) {}
124 124

	
125 125
      /// \brief Run the algorithm. It returns \c true if a path is found
126 126
      /// from the source node to the target node.
127 127
      bool run() {
128 128
        HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
129 129
        Heap heap(heap_cross_ref);
130 130
        heap.push(_s, 0);
131 131
        _pred[_s] = INVALID;
132 132
        _proc_nodes.clear();
133 133

	
134 134
        // Process nodes
135 135
        while (!heap.empty() && heap.top() != _t) {
136 136
          Node u = heap.top(), v;
137 137
          Length d = heap.prio() + _potential[u], nd;
138 138
          _dist[u] = heap.prio();
139 139
          heap.pop();
140 140
          _proc_nodes.push_back(u);
141 141

	
142 142
          // Traverse outgoing arcs
143 143
          for (OutArcIt e(_graph, u); e != INVALID; ++e) {
144 144
            if (_flow[e] == 0) {
145 145
              v = _graph.target(e);
146 146
              switch(heap.state(v)) {
147 147
              case Heap::PRE_HEAP:
148 148
                heap.push(v, d + _length[e] - _potential[v]);
149 149
                _pred[v] = e;
150 150
                break;
151 151
              case Heap::IN_HEAP:
152 152
                nd = d + _length[e] - _potential[v];
153 153
                if (nd < heap[v]) {
154 154
                  heap.decrease(v, nd);
155 155
                  _pred[v] = e;
156 156
                }
157 157
                break;
158 158
              case Heap::POST_HEAP:
159 159
                break;
160 160
              }
161 161
            }
162 162
          }
163 163

	
164 164
          // Traverse incoming arcs
165 165
          for (InArcIt e(_graph, u); e != INVALID; ++e) {
166 166
            if (_flow[e] == 1) {
167 167
              v = _graph.source(e);
168 168
              switch(heap.state(v)) {
169 169
              case Heap::PRE_HEAP:
170 170
                heap.push(v, d - _length[e] - _potential[v]);
171 171
                _pred[v] = e;
172 172
                break;
173 173
              case Heap::IN_HEAP:
174 174
                nd = d - _length[e] - _potential[v];
175 175
                if (nd < heap[v]) {
176 176
                  heap.decrease(v, nd);
177 177
                  _pred[v] = e;
178 178
                }
179 179
                break;
180 180
              case Heap::POST_HEAP:
181 181
                break;
182 182
              }
183 183
            }
184 184
          }
185 185
        }
186 186
        if (heap.empty()) return false;
187 187

	
188 188
        // Update potentials of processed nodes
189 189
        Length t_dist = heap.prio();
190 190
        for (int i = 0; i < int(_proc_nodes.size()); ++i)
191 191
          _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
192 192
        return true;
193 193
      }
194 194

	
195 195
    }; //class ResidualDijkstra
196 196

	
197 197
  private:
198 198

	
199 199
    // The digraph the algorithm runs on
200 200
    const Digraph &_graph;
201 201
    // The length map
202 202
    const LengthMap &_length;
203 203
    
204 204
    // Arc map of the current flow
205 205
    FlowMap *_flow;
206 206
    bool _local_flow;
207 207
    // Node map of the current potentials
208 208
    PotentialMap *_potential;
209 209
    bool _local_potential;
210 210

	
211 211
    // The source node
212 212
    Node _source;
213 213
    // The target node
214 214
    Node _target;
215 215

	
216 216
    // Container to store the found paths
217 217
    std::vector< SimplePath<Digraph> > paths;
218 218
    int _path_num;
219 219

	
220 220
    // The pred arc map
221 221
    PredMap _pred;
222 222
    // Implementation of the Dijkstra algorithm for finding augmenting
223 223
    // shortest paths in the residual network
224 224
    ResidualDijkstra *_dijkstra;
225 225

	
226 226
  public:
227 227

	
228 228
    /// \brief Constructor.
229 229
    ///
230 230
    /// Constructor.
231 231
    ///
232 232
    /// \param digraph The digraph the algorithm runs on.
233 233
    /// \param length The length (cost) values of the arcs.
234 234
    /// \param s The source node.
235 235
    /// \param t The target node.
236 236
    Suurballe( const Digraph &digraph,
237 237
               const LengthMap &length,
238 238
               Node s, Node t ) :
239 239
      _graph(digraph), _length(length), _flow(0), _local_flow(false),
240 240
      _potential(0), _local_potential(false), _source(s), _target(t),
241 241
      _pred(digraph) {}
242 242

	
243 243
    /// Destructor.
244 244
    ~Suurballe() {
245 245
      if (_local_flow) delete _flow;
246 246
      if (_local_potential) delete _potential;
247 247
      delete _dijkstra;
248 248
    }
249 249

	
250 250
    /// \brief Set the flow map.
251 251
    ///
252 252
    /// This function sets the flow map.
253 253
    ///
254 254
    /// The found flow contains only 0 and 1 values. It is the union of
255 255
    /// the found arc-disjoint paths.
256 256
    ///
257 257
    /// \return \c (*this)
258 258
    Suurballe& flowMap(FlowMap &map) {
259 259
      if (_local_flow) {
260 260
        delete _flow;
261 261
        _local_flow = false;
262 262
      }
263 263
      _flow = &map;
264 264
      return *this;
265 265
    }
266 266

	
267 267
    /// \brief Set the potential map.
268 268
    ///
269 269
    /// This function sets the potential map.
270 270
    ///
271 271
    /// The potentials provide the dual solution of the underlying 
272 272
    /// minimum cost flow problem.
273 273
    ///
274 274
    /// \return \c (*this)
275 275
    Suurballe& potentialMap(PotentialMap &map) {
276 276
      if (_local_potential) {
277 277
        delete _potential;
278 278
        _local_potential = false;
279 279
      }
280 280
      _potential = &map;
281 281
      return *this;
282 282
    }
283 283

	
284 284
    /// \name Execution control
285 285
    /// The simplest way to execute the algorithm is to call the run()
286 286
    /// function.
287 287
    /// \n
288 288
    /// If you only need the flow that is the union of the found
289 289
    /// arc-disjoint paths, you may call init() and findFlow().
290 290

	
291 291
    /// @{
292 292

	
293 293
    /// \brief Run the algorithm.
294 294
    ///
295 295
    /// This function runs the algorithm.
296 296
    ///
297 297
    /// \param k The number of paths to be found.
298 298
    ///
299 299
    /// \return \c k if there are at least \c k arc-disjoint paths from
300 300
    /// \c s to \c t in the digraph. Otherwise it returns the number of
301 301
    /// arc-disjoint paths found.
302 302
    ///
303 303
    /// \note Apart from the return value, <tt>s.run(k)</tt> is just a
304 304
    /// shortcut of the following code.
305 305
    /// \code
306 306
    ///   s.init();
307 307
    ///   s.findFlow(k);
308 308
    ///   s.findPaths();
309 309
    /// \endcode
310 310
    int run(int k = 2) {
311 311
      init();
312 312
      findFlow(k);
313 313
      findPaths();
314 314
      return _path_num;
315 315
    }
316 316

	
317 317
    /// \brief Initialize the algorithm.
318 318
    ///
319 319
    /// This function initializes the algorithm.
320 320
    void init() {
321 321
      // Initialize maps
322 322
      if (!_flow) {
323 323
        _flow = new FlowMap(_graph);
324 324
        _local_flow = true;
325 325
      }
326 326
      if (!_potential) {
327 327
        _potential = new PotentialMap(_graph);
328 328
        _local_potential = true;
329 329
      }
330 330
      for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
331 331
      for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
332 332

	
333 333
      _dijkstra = new ResidualDijkstra( _graph, *_flow, _length, 
334 334
                                        *_potential, _pred,
335 335
                                        _source, _target );
336 336
    }
337 337

	
338 338
    /// \brief Execute the successive shortest path algorithm to find
339 339
    /// an optimal flow.
340 340
    ///
341 341
    /// This function executes the successive shortest path algorithm to
342 342
    /// find a minimum cost flow, which is the union of \c k or less
343 343
    /// arc-disjoint paths.
344 344
    ///
345 345
    /// \return \c k if there are at least \c k arc-disjoint paths from
346 346
    /// \c s to \c t in the digraph. Otherwise it returns the number of
347 347
    /// arc-disjoint paths found.
348 348
    ///
349 349
    /// \pre \ref init() must be called before using this function.
350 350
    int findFlow(int k = 2) {
351 351
      // Find shortest paths
352 352
      _path_num = 0;
353 353
      while (_path_num < k) {
354 354
        // Run Dijkstra
355 355
        if (!_dijkstra->run()) break;
356 356
        ++_path_num;
357 357

	
358 358
        // Set the flow along the found shortest path
359 359
        Node u = _target;
360 360
        Arc e;
361 361
        while ((e = _pred[u]) != INVALID) {
362 362
          if (u == _graph.target(e)) {
363 363
            (*_flow)[e] = 1;
364 364
            u = _graph.source(e);
365 365
          } else {
366 366
            (*_flow)[e] = 0;
367 367
            u = _graph.target(e);
368 368
          }
369 369
        }
370 370
      }
371 371
      return _path_num;
372 372
    }
373 373
    
374 374
    /// \brief Compute the paths from the flow.
375 375
    ///
376 376
    /// This function computes the paths from the flow.
377 377
    ///
378 378
    /// \pre \ref init() and \ref findFlow() must be called before using
379 379
    /// this function.
380 380
    void findPaths() {
381 381
      // Create the residual flow map (the union of the paths not found
382 382
      // so far)
383 383
      FlowMap res_flow(_graph);
384 384
      for(ArcIt a(_graph); a != INVALID; ++a) res_flow[a] = (*_flow)[a];
385 385

	
386 386
      paths.clear();
387 387
      paths.resize(_path_num);
388 388
      for (int i = 0; i < _path_num; ++i) {
389 389
        Node n = _source;
390 390
        while (n != _target) {
391 391
          OutArcIt e(_graph, n);
392 392
          for ( ; res_flow[e] == 0; ++e) ;
393 393
          n = _graph.target(e);
394 394
          paths[i].addBack(e);
395 395
          res_flow[e] = 0;
396 396
        }
397 397
      }
398 398
    }
399 399

	
400 400
    /// @}
401 401

	
402 402
    /// \name Query Functions
403 403
    /// The results of the algorithm can be obtained using these
404 404
    /// functions.
405 405
    /// \n The algorithm should be executed before using them.
406 406

	
407 407
    /// @{
408 408

	
409 409
    /// \brief Return a const reference to the arc map storing the
410 410
    /// found flow.
411 411
    ///
412 412
    /// This function returns a const reference to the arc map storing
413 413
    /// the flow that is the union of the found arc-disjoint paths.
414 414
    ///
415 415
    /// \pre \ref run() or \ref findFlow() must be called before using
416 416
    /// this function.
417 417
    const FlowMap& flowMap() const {
418 418
      return *_flow;
419 419
    }
420 420

	
421 421
    /// \brief Return a const reference to the node map storing the
422 422
    /// found potentials (the dual solution).
423 423
    ///
424 424
    /// This function returns a const reference to the node map storing
425 425
    /// the found potentials that provide the dual solution of the
426 426
    /// underlying minimum cost flow problem.
427 427
    ///
428 428
    /// \pre \ref run() or \ref findFlow() must be called before using
429 429
    /// this function.
430 430
    const PotentialMap& potentialMap() const {
431 431
      return *_potential;
432 432
    }
433 433

	
434 434
    /// \brief Return the flow on the given arc.
435 435
    ///
436 436
    /// This function returns the flow on the given arc.
437 437
    /// It is \c 1 if the arc is involved in one of the found paths,
438 438
    /// otherwise it is \c 0.
439 439
    ///
440 440
    /// \pre \ref run() or \ref findFlow() must be called before using
441 441
    /// this function.
442 442
    int flow(const Arc& arc) const {
443 443
      return (*_flow)[arc];
444 444
    }
445 445

	
446 446
    /// \brief Return the potential of the given node.
447 447
    ///
448 448
    /// This function returns the potential of the given node.
449 449
    ///
450 450
    /// \pre \ref run() or \ref findFlow() must be called before using
451 451
    /// this function.
452 452
    Length potential(const Node& node) const {
453 453
      return (*_potential)[node];
454 454
    }
455 455

	
456 456
    /// \brief Return the total length (cost) of the found paths (flow).
457 457
    ///
458 458
    /// This function returns the total length (cost) of the found paths
459 459
    /// (flow). The complexity of the function is \f$ O(e) \f$.
460 460
    ///
461 461
    /// \pre \ref run() or \ref findFlow() must be called before using
462 462
    /// this function.
463 463
    Length totalLength() const {
464 464
      Length c = 0;
465 465
      for (ArcIt e(_graph); e != INVALID; ++e)
466 466
        c += (*_flow)[e] * _length[e];
467 467
      return c;
468 468
    }
469 469

	
470 470
    /// \brief Return the number of the found paths.
471 471
    ///
472 472
    /// This function returns the number of the found paths.
473 473
    ///
474 474
    /// \pre \ref run() or \ref findFlow() must be called before using
475 475
    /// this function.
476 476
    int pathNum() const {
477 477
      return _path_num;
478 478
    }
479 479

	
480 480
    /// \brief Return a const reference to the specified path.
481 481
    ///
482 482
    /// This function returns a const reference to the specified path.
483 483
    ///
484 484
    /// \param i The function returns the \c i-th path.
485 485
    /// \c i must be between \c 0 and <tt>%pathNum()-1</tt>.
486 486
    ///
487 487
    /// \pre \ref run() or \ref findPaths() must be called before using
488 488
    /// this function.
489 489
    Path path(int i) const {
490 490
      return paths[i];
491 491
    }
492 492

	
493 493
    /// @}
494 494

	
495 495
  }; //class Suurballe
496 496

	
497 497
  ///@}
498 498

	
499 499
} //namespace lemon
500 500

	
501 501
#endif //LEMON_SUURBALLE_H
0 comments (0 inline)