gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Add missing 'const' for query functions of algorithms
0 4 0
default
4 files changed with 10 insertions and 10 deletions:
↑ Collapse diff ↑
Ignore white space 6 line context
... ...
@@ -1232,521 +1232,521 @@
1232 1232
    void start(const Node&) {}
1233 1233
    void reach(const Node&) {}
1234 1234
    void process(const Node&) {}
1235 1235
    void discover(const Arc&) {}
1236 1236
    void examine(const Arc&) {}
1237 1237

	
1238 1238
    template <typename _Visitor>
1239 1239
    struct Constraints {
1240 1240
      void constraints() {
1241 1241
        Arc arc;
1242 1242
        Node node;
1243 1243
        visitor.start(node);
1244 1244
        visitor.reach(node);
1245 1245
        visitor.process(node);
1246 1246
        visitor.discover(arc);
1247 1247
        visitor.examine(arc);
1248 1248
      }
1249 1249
      _Visitor& visitor;
1250 1250
    };
1251 1251
  };
1252 1252
#endif
1253 1253

	
1254 1254
  /// \brief Default traits class of BfsVisit class.
1255 1255
  ///
1256 1256
  /// Default traits class of BfsVisit class.
1257 1257
  /// \tparam _Digraph The type of the digraph the algorithm runs on.
1258 1258
  template<class _Digraph>
1259 1259
  struct BfsVisitDefaultTraits {
1260 1260

	
1261 1261
    /// \brief The type of the digraph the algorithm runs on.
1262 1262
    typedef _Digraph Digraph;
1263 1263

	
1264 1264
    /// \brief The type of the map that indicates which nodes are reached.
1265 1265
    ///
1266 1266
    /// The type of the map that indicates which nodes are reached.
1267 1267
    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
1268 1268
    typedef typename Digraph::template NodeMap<bool> ReachedMap;
1269 1269

	
1270 1270
    /// \brief Instantiates a ReachedMap.
1271 1271
    ///
1272 1272
    /// This function instantiates a ReachedMap.
1273 1273
    /// \param digraph is the digraph, to which
1274 1274
    /// we would like to define the ReachedMap.
1275 1275
    static ReachedMap *createReachedMap(const Digraph &digraph) {
1276 1276
      return new ReachedMap(digraph);
1277 1277
    }
1278 1278

	
1279 1279
  };
1280 1280

	
1281 1281
  /// \ingroup search
1282 1282
  ///
1283 1283
  /// \brief %BFS algorithm class with visitor interface.
1284 1284
  ///
1285 1285
  /// This class provides an efficient implementation of the %BFS algorithm
1286 1286
  /// with visitor interface.
1287 1287
  ///
1288 1288
  /// The %BfsVisit class provides an alternative interface to the Bfs
1289 1289
  /// class. It works with callback mechanism, the BfsVisit object calls
1290 1290
  /// the member functions of the \c Visitor class on every BFS event.
1291 1291
  ///
1292 1292
  /// This interface of the BFS algorithm should be used in special cases
1293 1293
  /// when extra actions have to be performed in connection with certain
1294 1294
  /// events of the BFS algorithm. Otherwise consider to use Bfs or bfs()
1295 1295
  /// instead.
1296 1296
  ///
1297 1297
  /// \tparam _Digraph The type of the digraph the algorithm runs on.
1298 1298
  /// The default value is
1299 1299
  /// \ref ListDigraph. The value of _Digraph is not used directly by
1300 1300
  /// \ref BfsVisit, it is only passed to \ref BfsVisitDefaultTraits.
1301 1301
  /// \tparam _Visitor The Visitor type that is used by the algorithm.
1302 1302
  /// \ref BfsVisitor "BfsVisitor<_Digraph>" is an empty visitor, which
1303 1303
  /// does not observe the BFS events. If you want to observe the BFS
1304 1304
  /// events, you should implement your own visitor class.
1305 1305
  /// \tparam _Traits Traits class to set various data types used by the
1306 1306
  /// algorithm. The default traits class is
1307 1307
  /// \ref BfsVisitDefaultTraits "BfsVisitDefaultTraits<_Digraph>".
1308 1308
  /// See \ref BfsVisitDefaultTraits for the documentation of
1309 1309
  /// a BFS visit traits class.
1310 1310
#ifdef DOXYGEN
1311 1311
  template <typename _Digraph, typename _Visitor, typename _Traits>
1312 1312
#else
1313 1313
  template <typename _Digraph = ListDigraph,
1314 1314
            typename _Visitor = BfsVisitor<_Digraph>,
1315 1315
            typename _Traits = BfsVisitDefaultTraits<_Digraph> >
1316 1316
#endif
1317 1317
  class BfsVisit {
1318 1318
  public:
1319 1319

	
1320 1320
    ///The traits class.
1321 1321
    typedef _Traits Traits;
1322 1322

	
1323 1323
    ///The type of the digraph the algorithm runs on.
1324 1324
    typedef typename Traits::Digraph Digraph;
1325 1325

	
1326 1326
    ///The visitor type used by the algorithm.
1327 1327
    typedef _Visitor Visitor;
1328 1328

	
1329 1329
    ///The type of the map that indicates which nodes are reached.
1330 1330
    typedef typename Traits::ReachedMap ReachedMap;
1331 1331

	
1332 1332
  private:
1333 1333

	
1334 1334
    typedef typename Digraph::Node Node;
1335 1335
    typedef typename Digraph::NodeIt NodeIt;
1336 1336
    typedef typename Digraph::Arc Arc;
1337 1337
    typedef typename Digraph::OutArcIt OutArcIt;
1338 1338

	
1339 1339
    //Pointer to the underlying digraph.
1340 1340
    const Digraph *_digraph;
1341 1341
    //Pointer to the visitor object.
1342 1342
    Visitor *_visitor;
1343 1343
    //Pointer to the map of reached status of the nodes.
1344 1344
    ReachedMap *_reached;
1345 1345
    //Indicates if _reached is locally allocated (true) or not.
1346 1346
    bool local_reached;
1347 1347

	
1348 1348
    std::vector<typename Digraph::Node> _list;
1349 1349
    int _list_front, _list_back;
1350 1350

	
1351 1351
    //Creates the maps if necessary.
1352 1352
    void create_maps() {
1353 1353
      if(!_reached) {
1354 1354
        local_reached = true;
1355 1355
        _reached = Traits::createReachedMap(*_digraph);
1356 1356
      }
1357 1357
    }
1358 1358

	
1359 1359
  protected:
1360 1360

	
1361 1361
    BfsVisit() {}
1362 1362

	
1363 1363
  public:
1364 1364

	
1365 1365
    typedef BfsVisit Create;
1366 1366

	
1367 1367
    /// \name Named Template Parameters
1368 1368

	
1369 1369
    ///@{
1370 1370
    template <class T>
1371 1371
    struct SetReachedMapTraits : public Traits {
1372 1372
      typedef T ReachedMap;
1373 1373
      static ReachedMap *createReachedMap(const Digraph &digraph) {
1374 1374
        LEMON_ASSERT(false, "ReachedMap is not initialized");
1375 1375
        return 0; // ignore warnings
1376 1376
      }
1377 1377
    };
1378 1378
    /// \brief \ref named-templ-param "Named parameter" for setting
1379 1379
    /// ReachedMap type.
1380 1380
    ///
1381 1381
    /// \ref named-templ-param "Named parameter" for setting ReachedMap type.
1382 1382
    template <class T>
1383 1383
    struct SetReachedMap : public BfsVisit< Digraph, Visitor,
1384 1384
                                            SetReachedMapTraits<T> > {
1385 1385
      typedef BfsVisit< Digraph, Visitor, SetReachedMapTraits<T> > Create;
1386 1386
    };
1387 1387
    ///@}
1388 1388

	
1389 1389
  public:
1390 1390

	
1391 1391
    /// \brief Constructor.
1392 1392
    ///
1393 1393
    /// Constructor.
1394 1394
    ///
1395 1395
    /// \param digraph The digraph the algorithm runs on.
1396 1396
    /// \param visitor The visitor object of the algorithm.
1397 1397
    BfsVisit(const Digraph& digraph, Visitor& visitor)
1398 1398
      : _digraph(&digraph), _visitor(&visitor),
1399 1399
        _reached(0), local_reached(false) {}
1400 1400

	
1401 1401
    /// \brief Destructor.
1402 1402
    ~BfsVisit() {
1403 1403
      if(local_reached) delete _reached;
1404 1404
    }
1405 1405

	
1406 1406
    /// \brief Sets the map that indicates which nodes are reached.
1407 1407
    ///
1408 1408
    /// Sets the map that indicates which nodes are reached.
1409 1409
    /// If you don't use this function before calling \ref run(Node) "run()"
1410 1410
    /// or \ref init(), an instance will be allocated automatically.
1411 1411
    /// The destructor deallocates this automatically allocated map,
1412 1412
    /// of course.
1413 1413
    /// \return <tt> (*this) </tt>
1414 1414
    BfsVisit &reachedMap(ReachedMap &m) {
1415 1415
      if(local_reached) {
1416 1416
        delete _reached;
1417 1417
        local_reached = false;
1418 1418
      }
1419 1419
      _reached = &m;
1420 1420
      return *this;
1421 1421
    }
1422 1422

	
1423 1423
  public:
1424 1424

	
1425 1425
    /// \name Execution Control
1426 1426
    /// The simplest way to execute the BFS algorithm is to use one of the
1427 1427
    /// member functions called \ref run(Node) "run()".\n
1428 1428
    /// If you need more control on the execution, first you have to call
1429 1429
    /// \ref init(), then you can add several source nodes with
1430 1430
    /// \ref addSource(). Finally the actual path computation can be
1431 1431
    /// performed with one of the \ref start() functions.
1432 1432

	
1433 1433
    /// @{
1434 1434

	
1435 1435
    /// \brief Initializes the internal data structures.
1436 1436
    ///
1437 1437
    /// Initializes the internal data structures.
1438 1438
    void init() {
1439 1439
      create_maps();
1440 1440
      _list.resize(countNodes(*_digraph));
1441 1441
      _list_front = _list_back = -1;
1442 1442
      for (NodeIt u(*_digraph) ; u != INVALID ; ++u) {
1443 1443
        _reached->set(u, false);
1444 1444
      }
1445 1445
    }
1446 1446

	
1447 1447
    /// \brief Adds a new source node.
1448 1448
    ///
1449 1449
    /// Adds a new source node to the set of nodes to be processed.
1450 1450
    void addSource(Node s) {
1451 1451
      if(!(*_reached)[s]) {
1452 1452
          _reached->set(s,true);
1453 1453
          _visitor->start(s);
1454 1454
          _visitor->reach(s);
1455 1455
          _list[++_list_back] = s;
1456 1456
        }
1457 1457
    }
1458 1458

	
1459 1459
    /// \brief Processes the next node.
1460 1460
    ///
1461 1461
    /// Processes the next node.
1462 1462
    ///
1463 1463
    /// \return The processed node.
1464 1464
    ///
1465 1465
    /// \pre The queue must not be empty.
1466 1466
    Node processNextNode() {
1467 1467
      Node n = _list[++_list_front];
1468 1468
      _visitor->process(n);
1469 1469
      Arc e;
1470 1470
      for (_digraph->firstOut(e, n); e != INVALID; _digraph->nextOut(e)) {
1471 1471
        Node m = _digraph->target(e);
1472 1472
        if (!(*_reached)[m]) {
1473 1473
          _visitor->discover(e);
1474 1474
          _visitor->reach(m);
1475 1475
          _reached->set(m, true);
1476 1476
          _list[++_list_back] = m;
1477 1477
        } else {
1478 1478
          _visitor->examine(e);
1479 1479
        }
1480 1480
      }
1481 1481
      return n;
1482 1482
    }
1483 1483

	
1484 1484
    /// \brief Processes the next node.
1485 1485
    ///
1486 1486
    /// Processes the next node and checks if the given target node
1487 1487
    /// is reached. If the target node is reachable from the processed
1488 1488
    /// node, then the \c reach parameter will be set to \c true.
1489 1489
    ///
1490 1490
    /// \param target The target node.
1491 1491
    /// \retval reach Indicates if the target node is reached.
1492 1492
    /// It should be initially \c false.
1493 1493
    ///
1494 1494
    /// \return The processed node.
1495 1495
    ///
1496 1496
    /// \pre The queue must not be empty.
1497 1497
    Node processNextNode(Node target, bool& reach) {
1498 1498
      Node n = _list[++_list_front];
1499 1499
      _visitor->process(n);
1500 1500
      Arc e;
1501 1501
      for (_digraph->firstOut(e, n); e != INVALID; _digraph->nextOut(e)) {
1502 1502
        Node m = _digraph->target(e);
1503 1503
        if (!(*_reached)[m]) {
1504 1504
          _visitor->discover(e);
1505 1505
          _visitor->reach(m);
1506 1506
          _reached->set(m, true);
1507 1507
          _list[++_list_back] = m;
1508 1508
          reach = reach || (target == m);
1509 1509
        } else {
1510 1510
          _visitor->examine(e);
1511 1511
        }
1512 1512
      }
1513 1513
      return n;
1514 1514
    }
1515 1515

	
1516 1516
    /// \brief Processes the next node.
1517 1517
    ///
1518 1518
    /// Processes the next node and checks if at least one of reached
1519 1519
    /// nodes has \c true value in the \c nm node map. If one node
1520 1520
    /// with \c true value is reachable from the processed node, then the
1521 1521
    /// \c rnode parameter will be set to the first of such nodes.
1522 1522
    ///
1523 1523
    /// \param nm A \c bool (or convertible) node map that indicates the
1524 1524
    /// possible targets.
1525 1525
    /// \retval rnode The reached target node.
1526 1526
    /// It should be initially \c INVALID.
1527 1527
    ///
1528 1528
    /// \return The processed node.
1529 1529
    ///
1530 1530
    /// \pre The queue must not be empty.
1531 1531
    template <typename NM>
1532 1532
    Node processNextNode(const NM& nm, Node& rnode) {
1533 1533
      Node n = _list[++_list_front];
1534 1534
      _visitor->process(n);
1535 1535
      Arc e;
1536 1536
      for (_digraph->firstOut(e, n); e != INVALID; _digraph->nextOut(e)) {
1537 1537
        Node m = _digraph->target(e);
1538 1538
        if (!(*_reached)[m]) {
1539 1539
          _visitor->discover(e);
1540 1540
          _visitor->reach(m);
1541 1541
          _reached->set(m, true);
1542 1542
          _list[++_list_back] = m;
1543 1543
          if (nm[m] && rnode == INVALID) rnode = m;
1544 1544
        } else {
1545 1545
          _visitor->examine(e);
1546 1546
        }
1547 1547
      }
1548 1548
      return n;
1549 1549
    }
1550 1550

	
1551 1551
    /// \brief The next node to be processed.
1552 1552
    ///
1553 1553
    /// Returns the next node to be processed or \c INVALID if the queue
1554 1554
    /// is empty.
1555 1555
    Node nextNode() const {
1556 1556
      return _list_front != _list_back ? _list[_list_front + 1] : INVALID;
1557 1557
    }
1558 1558

	
1559 1559
    /// \brief Returns \c false if there are nodes
1560 1560
    /// to be processed.
1561 1561
    ///
1562 1562
    /// Returns \c false if there are nodes
1563 1563
    /// to be processed in the queue.
1564 1564
    bool emptyQueue() const { return _list_front == _list_back; }
1565 1565

	
1566 1566
    /// \brief Returns the number of the nodes to be processed.
1567 1567
    ///
1568 1568
    /// Returns the number of the nodes to be processed in the queue.
1569 1569
    int queueSize() const { return _list_back - _list_front; }
1570 1570

	
1571 1571
    /// \brief Executes the algorithm.
1572 1572
    ///
1573 1573
    /// Executes the algorithm.
1574 1574
    ///
1575 1575
    /// This method runs the %BFS algorithm from the root node(s)
1576 1576
    /// in order to compute the shortest path to each node.
1577 1577
    ///
1578 1578
    /// The algorithm computes
1579 1579
    /// - the shortest path tree (forest),
1580 1580
    /// - the distance of each node from the root(s).
1581 1581
    ///
1582 1582
    /// \pre init() must be called and at least one root node should be added
1583 1583
    /// with addSource() before using this function.
1584 1584
    ///
1585 1585
    /// \note <tt>b.start()</tt> is just a shortcut of the following code.
1586 1586
    /// \code
1587 1587
    ///   while ( !b.emptyQueue() ) {
1588 1588
    ///     b.processNextNode();
1589 1589
    ///   }
1590 1590
    /// \endcode
1591 1591
    void start() {
1592 1592
      while ( !emptyQueue() ) processNextNode();
1593 1593
    }
1594 1594

	
1595 1595
    /// \brief Executes the algorithm until the given target node is reached.
1596 1596
    ///
1597 1597
    /// Executes the algorithm until the given target node is reached.
1598 1598
    ///
1599 1599
    /// This method runs the %BFS algorithm from the root node(s)
1600 1600
    /// in order to compute the shortest path to \c t.
1601 1601
    ///
1602 1602
    /// The algorithm computes
1603 1603
    /// - the shortest path to \c t,
1604 1604
    /// - the distance of \c t from the root(s).
1605 1605
    ///
1606 1606
    /// \pre init() must be called and at least one root node should be
1607 1607
    /// added with addSource() before using this function.
1608 1608
    ///
1609 1609
    /// \note <tt>b.start(t)</tt> is just a shortcut of the following code.
1610 1610
    /// \code
1611 1611
    ///   bool reach = false;
1612 1612
    ///   while ( !b.emptyQueue() && !reach ) {
1613 1613
    ///     b.processNextNode(t, reach);
1614 1614
    ///   }
1615 1615
    /// \endcode
1616 1616
    void start(Node t) {
1617 1617
      bool reach = false;
1618 1618
      while ( !emptyQueue() && !reach ) processNextNode(t, reach);
1619 1619
    }
1620 1620

	
1621 1621
    /// \brief Executes the algorithm until a condition is met.
1622 1622
    ///
1623 1623
    /// Executes the algorithm until a condition is met.
1624 1624
    ///
1625 1625
    /// This method runs the %BFS algorithm from the root node(s) in
1626 1626
    /// order to compute the shortest path to a node \c v with
1627 1627
    /// <tt>nm[v]</tt> true, if such a node can be found.
1628 1628
    ///
1629 1629
    /// \param nm must be a bool (or convertible) node map. The
1630 1630
    /// algorithm will stop when it reaches a node \c v with
1631 1631
    /// <tt>nm[v]</tt> true.
1632 1632
    ///
1633 1633
    /// \return The reached node \c v with <tt>nm[v]</tt> true or
1634 1634
    /// \c INVALID if no such node was found.
1635 1635
    ///
1636 1636
    /// \pre init() must be called and at least one root node should be
1637 1637
    /// added with addSource() before using this function.
1638 1638
    ///
1639 1639
    /// \note <tt>b.start(nm)</tt> is just a shortcut of the following code.
1640 1640
    /// \code
1641 1641
    ///   Node rnode = INVALID;
1642 1642
    ///   while ( !b.emptyQueue() && rnode == INVALID ) {
1643 1643
    ///     b.processNextNode(nm, rnode);
1644 1644
    ///   }
1645 1645
    ///   return rnode;
1646 1646
    /// \endcode
1647 1647
    template <typename NM>
1648 1648
    Node start(const NM &nm) {
1649 1649
      Node rnode = INVALID;
1650 1650
      while ( !emptyQueue() && rnode == INVALID ) {
1651 1651
        processNextNode(nm, rnode);
1652 1652
      }
1653 1653
      return rnode;
1654 1654
    }
1655 1655

	
1656 1656
    /// \brief Runs the algorithm from the given source node.
1657 1657
    ///
1658 1658
    /// This method runs the %BFS algorithm from node \c s
1659 1659
    /// in order to compute the shortest path to each node.
1660 1660
    ///
1661 1661
    /// The algorithm computes
1662 1662
    /// - the shortest path tree,
1663 1663
    /// - the distance of each node from the root.
1664 1664
    ///
1665 1665
    /// \note <tt>b.run(s)</tt> is just a shortcut of the following code.
1666 1666
    ///\code
1667 1667
    ///   b.init();
1668 1668
    ///   b.addSource(s);
1669 1669
    ///   b.start();
1670 1670
    ///\endcode
1671 1671
    void run(Node s) {
1672 1672
      init();
1673 1673
      addSource(s);
1674 1674
      start();
1675 1675
    }
1676 1676

	
1677 1677
    /// \brief Finds the shortest path between \c s and \c t.
1678 1678
    ///
1679 1679
    /// This method runs the %BFS algorithm from node \c s
1680 1680
    /// in order to compute the shortest path to node \c t
1681 1681
    /// (it stops searching when \c t is processed).
1682 1682
    ///
1683 1683
    /// \return \c true if \c t is reachable form \c s.
1684 1684
    ///
1685 1685
    /// \note Apart from the return value, <tt>b.run(s,t)</tt> is just a
1686 1686
    /// shortcut of the following code.
1687 1687
    ///\code
1688 1688
    ///   b.init();
1689 1689
    ///   b.addSource(s);
1690 1690
    ///   b.start(t);
1691 1691
    ///\endcode
1692 1692
    bool run(Node s,Node t) {
1693 1693
      init();
1694 1694
      addSource(s);
1695 1695
      start(t);
1696 1696
      return reached(t);
1697 1697
    }
1698 1698

	
1699 1699
    /// \brief Runs the algorithm to visit all nodes in the digraph.
1700 1700
    ///
1701 1701
    /// This method runs the %BFS algorithm in order to
1702 1702
    /// compute the shortest path to each node.
1703 1703
    ///
1704 1704
    /// The algorithm computes
1705 1705
    /// - the shortest path tree (forest),
1706 1706
    /// - the distance of each node from the root(s).
1707 1707
    ///
1708 1708
    /// \note <tt>b.run(s)</tt> is just a shortcut of the following code.
1709 1709
    ///\code
1710 1710
    ///  b.init();
1711 1711
    ///  for (NodeIt n(gr); n != INVALID; ++n) {
1712 1712
    ///    if (!b.reached(n)) {
1713 1713
    ///      b.addSource(n);
1714 1714
    ///      b.start();
1715 1715
    ///    }
1716 1716
    ///  }
1717 1717
    ///\endcode
1718 1718
    void run() {
1719 1719
      init();
1720 1720
      for (NodeIt it(*_digraph); it != INVALID; ++it) {
1721 1721
        if (!reached(it)) {
1722 1722
          addSource(it);
1723 1723
          start();
1724 1724
        }
1725 1725
      }
1726 1726
    }
1727 1727

	
1728 1728
    ///@}
1729 1729

	
1730 1730
    /// \name Query Functions
1731 1731
    /// The results of the BFS algorithm can be obtained using these
1732 1732
    /// functions.\n
1733 1733
    /// Either \ref run(Node) "run()" or \ref start() should be called
1734 1734
    /// before using them.
1735 1735

	
1736 1736
    ///@{
1737 1737

	
1738 1738
    /// \brief Checks if a node is reached from the root(s).
1739 1739
    ///
1740 1740
    /// Returns \c true if \c v is reached from the root(s).
1741 1741
    ///
1742 1742
    /// \pre Either \ref run(Node) "run()" or \ref init()
1743 1743
    /// must be called before using this function.
1744
    bool reached(Node v) { return (*_reached)[v]; }
1744
    bool reached(Node v) const { return (*_reached)[v]; }
1745 1745

	
1746 1746
    ///@}
1747 1747

	
1748 1748
  };
1749 1749

	
1750 1750
} //END OF NAMESPACE LEMON
1751 1751

	
1752 1752
#endif
Ignore white space 1024 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_CIRCULATION_H
20 20
#define LEMON_CIRCULATION_H
21 21

	
22 22
#include <lemon/tolerance.h>
23 23
#include <lemon/elevator.h>
24 24

	
25 25
///\ingroup max_flow
26 26
///\file
27 27
///\brief Push-relabel algorithm for finding a feasible circulation.
28 28
///
29 29
namespace lemon {
30 30

	
31 31
  /// \brief Default traits class of Circulation class.
32 32
  ///
33 33
  /// Default traits class of Circulation class.
34 34
  /// \tparam _Diraph Digraph type.
35 35
  /// \tparam _LCapMap Lower bound capacity map type.
36 36
  /// \tparam _UCapMap Upper bound capacity map type.
37 37
  /// \tparam _DeltaMap Delta map type.
38 38
  template <typename _Diraph, typename _LCapMap,
39 39
            typename _UCapMap, typename _DeltaMap>
40 40
  struct CirculationDefaultTraits {
41 41

	
42 42
    /// \brief The type of the digraph the algorithm runs on.
43 43
    typedef _Diraph Digraph;
44 44

	
45 45
    /// \brief The type of the map that stores the circulation lower
46 46
    /// bound.
47 47
    ///
48 48
    /// The type of the map that stores the circulation lower bound.
49 49
    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
50 50
    typedef _LCapMap LCapMap;
51 51

	
52 52
    /// \brief The type of the map that stores the circulation upper
53 53
    /// bound.
54 54
    ///
55 55
    /// The type of the map that stores the circulation upper bound.
56 56
    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
57 57
    typedef _UCapMap UCapMap;
58 58

	
59 59
    /// \brief The type of the map that stores the lower bound for
60 60
    /// the supply of the nodes.
61 61
    ///
62 62
    /// The type of the map that stores the lower bound for the supply
63 63
    /// of the nodes. It must meet the \ref concepts::ReadMap "ReadMap"
64 64
    /// concept.
65 65
    typedef _DeltaMap DeltaMap;
66 66

	
67 67
    /// \brief The type of the flow values.
68 68
    typedef typename DeltaMap::Value Value;
69 69

	
70 70
    /// \brief The type of the map that stores the flow values.
71 71
    ///
72 72
    /// The type of the map that stores the flow values.
73 73
    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
74 74
    typedef typename Digraph::template ArcMap<Value> FlowMap;
75 75

	
76 76
    /// \brief Instantiates a FlowMap.
77 77
    ///
78 78
    /// This function instantiates a \ref FlowMap.
79 79
    /// \param digraph The digraph, to which we would like to define
80 80
    /// the flow map.
81 81
    static FlowMap* createFlowMap(const Digraph& digraph) {
82 82
      return new FlowMap(digraph);
83 83
    }
84 84

	
85 85
    /// \brief The elevator type used by the algorithm.
86 86
    ///
87 87
    /// The elevator type used by the algorithm.
88 88
    ///
89 89
    /// \sa Elevator
90 90
    /// \sa LinkedElevator
91 91
    typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
92 92

	
93 93
    /// \brief Instantiates an Elevator.
94 94
    ///
95 95
    /// This function instantiates an \ref Elevator.
96 96
    /// \param digraph The digraph, to which we would like to define
97 97
    /// the elevator.
98 98
    /// \param max_level The maximum level of the elevator.
99 99
    static Elevator* createElevator(const Digraph& digraph, int max_level) {
100 100
      return new Elevator(digraph, max_level);
101 101
    }
102 102

	
103 103
    /// \brief The tolerance used by the algorithm
104 104
    ///
105 105
    /// The tolerance used by the algorithm to handle inexact computation.
106 106
    typedef lemon::Tolerance<Value> Tolerance;
107 107

	
108 108
  };
109 109

	
110 110
  /**
111 111
     \brief Push-relabel algorithm for the network circulation problem.
112 112

	
113 113
     \ingroup max_flow
114 114
     This class implements a push-relabel algorithm for the network
115 115
     circulation problem.
116 116
     It is to find a feasible circulation when lower and upper bounds
117 117
     are given for the flow values on the arcs and lower bounds
118 118
     are given for the supply values of the nodes.
119 119

	
120 120
     The exact formulation of this problem is the following.
121 121
     Let \f$G=(V,A)\f$ be a digraph,
122 122
     \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$,
123 123
     \f$delta: V\rightarrow\mathbf{R}\f$. Find a feasible circulation
124 124
     \f$f: A\rightarrow\mathbf{R}^+_0\f$ so that
125 125
     \f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a)
126 126
     \geq delta(v) \quad \forall v\in V, \f]
127 127
     \f[ lower(a)\leq f(a) \leq upper(a) \quad \forall a\in A. \f]
128 128
     \note \f$delta(v)\f$ specifies a lower bound for the supply of node
129 129
     \f$v\f$. It can be either positive or negative, however note that
130 130
     \f$\sum_{v\in V}delta(v)\f$ should be zero or negative in order to
131 131
     have a feasible solution.
132 132

	
133 133
     \note A special case of this problem is when
134 134
     \f$\sum_{v\in V}delta(v) = 0\f$. Then the supply of each node \f$v\f$
135 135
     will be \e equal \e to \f$delta(v)\f$, if a circulation can be found.
136 136
     Thus a feasible solution for the
137 137
     \ref min_cost_flow "minimum cost flow" problem can be calculated
138 138
     in this way.
139 139

	
140 140
     \tparam _Digraph The type of the digraph the algorithm runs on.
141 141
     \tparam _LCapMap The type of the lower bound capacity map. The default
142 142
     map type is \ref concepts::Digraph::ArcMap "_Digraph::ArcMap<int>".
143 143
     \tparam _UCapMap The type of the upper bound capacity map. The default
144 144
     map type is \c _LCapMap.
145 145
     \tparam _DeltaMap The type of the map that stores the lower bound
146 146
     for the supply of the nodes. The default map type is
147 147
     \c _Digraph::ArcMap<_UCapMap::Value>.
148 148
  */
149 149
#ifdef DOXYGEN
150 150
template< typename _Digraph,
151 151
          typename _LCapMap,
152 152
          typename _UCapMap,
153 153
          typename _DeltaMap,
154 154
          typename _Traits >
155 155
#else
156 156
template< typename _Digraph,
157 157
          typename _LCapMap = typename _Digraph::template ArcMap<int>,
158 158
          typename _UCapMap = _LCapMap,
159 159
          typename _DeltaMap = typename _Digraph::
160 160
                               template NodeMap<typename _UCapMap::Value>,
161 161
          typename _Traits=CirculationDefaultTraits<_Digraph, _LCapMap,
162 162
                                                    _UCapMap, _DeltaMap> >
163 163
#endif
164 164
  class Circulation {
165 165
  public:
166 166

	
167 167
    ///The \ref CirculationDefaultTraits "traits class" of the algorithm.
168 168
    typedef _Traits Traits;
169 169
    ///The type of the digraph the algorithm runs on.
170 170
    typedef typename Traits::Digraph Digraph;
171 171
    ///The type of the flow values.
172 172
    typedef typename Traits::Value Value;
173 173

	
174 174
    /// The type of the lower bound capacity map.
175 175
    typedef typename Traits::LCapMap LCapMap;
176 176
    /// The type of the upper bound capacity map.
177 177
    typedef typename Traits::UCapMap UCapMap;
178 178
    /// \brief The type of the map that stores the lower bound for
179 179
    /// the supply of the nodes.
180 180
    typedef typename Traits::DeltaMap DeltaMap;
181 181
    ///The type of the flow map.
182 182
    typedef typename Traits::FlowMap FlowMap;
183 183

	
184 184
    ///The type of the elevator.
185 185
    typedef typename Traits::Elevator Elevator;
186 186
    ///The type of the tolerance.
187 187
    typedef typename Traits::Tolerance Tolerance;
188 188

	
189 189
  private:
190 190

	
191 191
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
192 192

	
193 193
    const Digraph &_g;
194 194
    int _node_num;
195 195

	
196 196
    const LCapMap *_lo;
197 197
    const UCapMap *_up;
198 198
    const DeltaMap *_delta;
199 199

	
200 200
    FlowMap *_flow;
201 201
    bool _local_flow;
202 202

	
203 203
    Elevator* _level;
204 204
    bool _local_level;
205 205

	
206 206
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
207 207
    ExcessMap* _excess;
208 208

	
209 209
    Tolerance _tol;
210 210
    int _el;
211 211

	
212 212
  public:
213 213

	
214 214
    typedef Circulation Create;
215 215

	
216 216
    ///\name Named Template Parameters
217 217

	
218 218
    ///@{
219 219

	
220 220
    template <typename _FlowMap>
221 221
    struct SetFlowMapTraits : public Traits {
222 222
      typedef _FlowMap FlowMap;
223 223
      static FlowMap *createFlowMap(const Digraph&) {
224 224
        LEMON_ASSERT(false, "FlowMap is not initialized");
225 225
        return 0; // ignore warnings
226 226
      }
227 227
    };
228 228

	
229 229
    /// \brief \ref named-templ-param "Named parameter" for setting
230 230
    /// FlowMap type
231 231
    ///
232 232
    /// \ref named-templ-param "Named parameter" for setting FlowMap
233 233
    /// type.
234 234
    template <typename _FlowMap>
235 235
    struct SetFlowMap
236 236
      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
237 237
                           SetFlowMapTraits<_FlowMap> > {
238 238
      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
239 239
                          SetFlowMapTraits<_FlowMap> > Create;
240 240
    };
241 241

	
242 242
    template <typename _Elevator>
243 243
    struct SetElevatorTraits : public Traits {
244 244
      typedef _Elevator Elevator;
245 245
      static Elevator *createElevator(const Digraph&, int) {
246 246
        LEMON_ASSERT(false, "Elevator is not initialized");
247 247
        return 0; // ignore warnings
248 248
      }
249 249
    };
250 250

	
251 251
    /// \brief \ref named-templ-param "Named parameter" for setting
252 252
    /// Elevator type
253 253
    ///
254 254
    /// \ref named-templ-param "Named parameter" for setting Elevator
255 255
    /// type. If this named parameter is used, then an external
256 256
    /// elevator object must be passed to the algorithm using the
257 257
    /// \ref elevator(Elevator&) "elevator()" function before calling
258 258
    /// \ref run() or \ref init().
259 259
    /// \sa SetStandardElevator
260 260
    template <typename _Elevator>
261 261
    struct SetElevator
262 262
      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
263 263
                           SetElevatorTraits<_Elevator> > {
264 264
      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
265 265
                          SetElevatorTraits<_Elevator> > Create;
266 266
    };
267 267

	
268 268
    template <typename _Elevator>
269 269
    struct SetStandardElevatorTraits : public Traits {
270 270
      typedef _Elevator Elevator;
271 271
      static Elevator *createElevator(const Digraph& digraph, int max_level) {
272 272
        return new Elevator(digraph, max_level);
273 273
      }
274 274
    };
275 275

	
276 276
    /// \brief \ref named-templ-param "Named parameter" for setting
277 277
    /// Elevator type with automatic allocation
278 278
    ///
279 279
    /// \ref named-templ-param "Named parameter" for setting Elevator
280 280
    /// type with automatic allocation.
281 281
    /// The Elevator should have standard constructor interface to be
282 282
    /// able to automatically created by the algorithm (i.e. the
283 283
    /// digraph and the maximum level should be passed to it).
284 284
    /// However an external elevator object could also be passed to the
285 285
    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
286 286
    /// before calling \ref run() or \ref init().
287 287
    /// \sa SetElevator
288 288
    template <typename _Elevator>
289 289
    struct SetStandardElevator
290 290
      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
291 291
                       SetStandardElevatorTraits<_Elevator> > {
292 292
      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
293 293
                      SetStandardElevatorTraits<_Elevator> > Create;
294 294
    };
295 295

	
296 296
    /// @}
297 297

	
298 298
  protected:
299 299

	
300 300
    Circulation() {}
301 301

	
302 302
  public:
303 303

	
304 304
    /// The constructor of the class.
305 305

	
306 306
    /// The constructor of the class.
307 307
    /// \param g The digraph the algorithm runs on.
308 308
    /// \param lo The lower bound capacity of the arcs.
309 309
    /// \param up The upper bound capacity of the arcs.
310 310
    /// \param delta The lower bound for the supply of the nodes.
311 311
    Circulation(const Digraph &g,const LCapMap &lo,
312 312
                const UCapMap &up,const DeltaMap &delta)
313 313
      : _g(g), _node_num(),
314 314
        _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
315 315
        _level(0), _local_level(false), _excess(0), _el() {}
316 316

	
317 317
    /// Destructor.
318 318
    ~Circulation() {
319 319
      destroyStructures();
320 320
    }
321 321

	
322 322

	
323 323
  private:
324 324

	
325 325
    void createStructures() {
326 326
      _node_num = _el = countNodes(_g);
327 327

	
328 328
      if (!_flow) {
329 329
        _flow = Traits::createFlowMap(_g);
330 330
        _local_flow = true;
331 331
      }
332 332
      if (!_level) {
333 333
        _level = Traits::createElevator(_g, _node_num);
334 334
        _local_level = true;
335 335
      }
336 336
      if (!_excess) {
337 337
        _excess = new ExcessMap(_g);
338 338
      }
339 339
    }
340 340

	
341 341
    void destroyStructures() {
342 342
      if (_local_flow) {
343 343
        delete _flow;
344 344
      }
345 345
      if (_local_level) {
346 346
        delete _level;
347 347
      }
348 348
      if (_excess) {
349 349
        delete _excess;
350 350
      }
351 351
    }
352 352

	
353 353
  public:
354 354

	
355 355
    /// Sets the lower bound capacity map.
356 356

	
357 357
    /// Sets the lower bound capacity map.
358 358
    /// \return <tt>(*this)</tt>
359 359
    Circulation& lowerCapMap(const LCapMap& map) {
360 360
      _lo = &map;
361 361
      return *this;
362 362
    }
363 363

	
364 364
    /// Sets the upper bound capacity map.
365 365

	
366 366
    /// Sets the upper bound capacity map.
367 367
    /// \return <tt>(*this)</tt>
368 368
    Circulation& upperCapMap(const LCapMap& map) {
369 369
      _up = &map;
370 370
      return *this;
371 371
    }
372 372

	
373 373
    /// Sets the lower bound map for the supply of the nodes.
374 374

	
375 375
    /// Sets the lower bound map for the supply of the nodes.
376 376
    /// \return <tt>(*this)</tt>
377 377
    Circulation& deltaMap(const DeltaMap& map) {
378 378
      _delta = &map;
379 379
      return *this;
380 380
    }
381 381

	
382 382
    /// \brief Sets the flow map.
383 383
    ///
384 384
    /// Sets the flow map.
385 385
    /// If you don't use this function before calling \ref run() or
386 386
    /// \ref init(), an instance will be allocated automatically.
387 387
    /// The destructor deallocates this automatically allocated map,
388 388
    /// of course.
389 389
    /// \return <tt>(*this)</tt>
390 390
    Circulation& flowMap(FlowMap& map) {
391 391
      if (_local_flow) {
392 392
        delete _flow;
393 393
        _local_flow = false;
394 394
      }
395 395
      _flow = &map;
396 396
      return *this;
397 397
    }
398 398

	
399 399
    /// \brief Sets the elevator used by algorithm.
400 400
    ///
401 401
    /// Sets the elevator used by algorithm.
402 402
    /// If you don't use this function before calling \ref run() or
403 403
    /// \ref init(), an instance will be allocated automatically.
404 404
    /// The destructor deallocates this automatically allocated elevator,
405 405
    /// of course.
406 406
    /// \return <tt>(*this)</tt>
407 407
    Circulation& elevator(Elevator& elevator) {
408 408
      if (_local_level) {
409 409
        delete _level;
410 410
        _local_level = false;
411 411
      }
412 412
      _level = &elevator;
413 413
      return *this;
414 414
    }
415 415

	
416 416
    /// \brief Returns a const reference to the elevator.
417 417
    ///
418 418
    /// Returns a const reference to the elevator.
419 419
    ///
420 420
    /// \pre Either \ref run() or \ref init() must be called before
421 421
    /// using this function.
422
    const Elevator& elevator() {
422
    const Elevator& elevator() const {
423 423
      return *_level;
424 424
    }
425 425

	
426 426
    /// \brief Sets the tolerance used by algorithm.
427 427
    ///
428 428
    /// Sets the tolerance used by algorithm.
429 429
    Circulation& tolerance(const Tolerance& tolerance) const {
430 430
      _tol = tolerance;
431 431
      return *this;
432 432
    }
433 433

	
434 434
    /// \brief Returns a const reference to the tolerance.
435 435
    ///
436 436
    /// Returns a const reference to the tolerance.
437 437
    const Tolerance& tolerance() const {
438 438
      return tolerance;
439 439
    }
440 440

	
441 441
    /// \name Execution Control
442 442
    /// The simplest way to execute the algorithm is to call \ref run().\n
443 443
    /// If you need more control on the initial solution or the execution,
444 444
    /// first you have to call one of the \ref init() functions, then
445 445
    /// the \ref start() function.
446 446

	
447 447
    ///@{
448 448

	
449 449
    /// Initializes the internal data structures.
450 450

	
451 451
    /// Initializes the internal data structures and sets all flow values
452 452
    /// to the lower bound.
453 453
    void init()
454 454
    {
455 455
      createStructures();
456 456

	
457 457
      for(NodeIt n(_g);n!=INVALID;++n) {
458 458
        _excess->set(n, (*_delta)[n]);
459 459
      }
460 460

	
461 461
      for (ArcIt e(_g);e!=INVALID;++e) {
462 462
        _flow->set(e, (*_lo)[e]);
463 463
        _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
464 464
        _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
465 465
      }
466 466

	
467 467
      // global relabeling tested, but in general case it provides
468 468
      // worse performance for random digraphs
469 469
      _level->initStart();
470 470
      for(NodeIt n(_g);n!=INVALID;++n)
471 471
        _level->initAddItem(n);
472 472
      _level->initFinish();
473 473
      for(NodeIt n(_g);n!=INVALID;++n)
474 474
        if(_tol.positive((*_excess)[n]))
475 475
          _level->activate(n);
476 476
    }
477 477

	
478 478
    /// Initializes the internal data structures using a greedy approach.
479 479

	
480 480
    /// Initializes the internal data structures using a greedy approach
481 481
    /// to construct the initial solution.
482 482
    void greedyInit()
483 483
    {
484 484
      createStructures();
485 485

	
486 486
      for(NodeIt n(_g);n!=INVALID;++n) {
487 487
        _excess->set(n, (*_delta)[n]);
488 488
      }
489 489

	
490 490
      for (ArcIt e(_g);e!=INVALID;++e) {
491 491
        if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
492 492
          _flow->set(e, (*_up)[e]);
493 493
          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
494 494
          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
495 495
        } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
496 496
          _flow->set(e, (*_lo)[e]);
497 497
          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_lo)[e]);
498 498
          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_lo)[e]);
499 499
        } else {
500 500
          Value fc = -(*_excess)[_g.target(e)];
501 501
          _flow->set(e, fc);
502 502
          _excess->set(_g.target(e), 0);
503 503
          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
504 504
        }
505 505
      }
506 506

	
507 507
      _level->initStart();
508 508
      for(NodeIt n(_g);n!=INVALID;++n)
509 509
        _level->initAddItem(n);
510 510
      _level->initFinish();
511 511
      for(NodeIt n(_g);n!=INVALID;++n)
512 512
        if(_tol.positive((*_excess)[n]))
513 513
          _level->activate(n);
514 514
    }
515 515

	
516 516
    ///Executes the algorithm
517 517

	
518 518
    ///This function executes the algorithm.
519 519
    ///
520 520
    ///\return \c true if a feasible circulation is found.
521 521
    ///
522 522
    ///\sa barrier()
523 523
    ///\sa barrierMap()
524 524
    bool start()
525 525
    {
526 526

	
527 527
      Node act;
528 528
      Node bact=INVALID;
529 529
      Node last_activated=INVALID;
530 530
      while((act=_level->highestActive())!=INVALID) {
531 531
        int actlevel=(*_level)[act];
532 532
        int mlevel=_node_num;
533 533
        Value exc=(*_excess)[act];
534 534

	
535 535
        for(OutArcIt e(_g,act);e!=INVALID; ++e) {
536 536
          Node v = _g.target(e);
537 537
          Value fc=(*_up)[e]-(*_flow)[e];
538 538
          if(!_tol.positive(fc)) continue;
539 539
          if((*_level)[v]<actlevel) {
540 540
            if(!_tol.less(fc, exc)) {
541 541
              _flow->set(e, (*_flow)[e] + exc);
542 542
              _excess->set(v, (*_excess)[v] + exc);
543 543
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
544 544
                _level->activate(v);
545 545
              _excess->set(act,0);
546 546
              _level->deactivate(act);
547 547
              goto next_l;
548 548
            }
549 549
            else {
550 550
              _flow->set(e, (*_up)[e]);
551 551
              _excess->set(v, (*_excess)[v] + fc);
552 552
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
553 553
                _level->activate(v);
554 554
              exc-=fc;
555 555
            }
556 556
          }
557 557
          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
558 558
        }
559 559
        for(InArcIt e(_g,act);e!=INVALID; ++e) {
560 560
          Node v = _g.source(e);
561 561
          Value fc=(*_flow)[e]-(*_lo)[e];
562 562
          if(!_tol.positive(fc)) continue;
563 563
          if((*_level)[v]<actlevel) {
564 564
            if(!_tol.less(fc, exc)) {
565 565
              _flow->set(e, (*_flow)[e] - exc);
566 566
              _excess->set(v, (*_excess)[v] + exc);
567 567
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
568 568
                _level->activate(v);
569 569
              _excess->set(act,0);
570 570
              _level->deactivate(act);
571 571
              goto next_l;
572 572
            }
573 573
            else {
574 574
              _flow->set(e, (*_lo)[e]);
575 575
              _excess->set(v, (*_excess)[v] + fc);
576 576
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
577 577
                _level->activate(v);
578 578
              exc-=fc;
579 579
            }
580 580
          }
581 581
          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
582 582
        }
583 583

	
584 584
        _excess->set(act, exc);
585 585
        if(!_tol.positive(exc)) _level->deactivate(act);
586 586
        else if(mlevel==_node_num) {
587 587
          _level->liftHighestActiveToTop();
588 588
          _el = _node_num;
589 589
          return false;
590 590
        }
591 591
        else {
592 592
          _level->liftHighestActive(mlevel+1);
593 593
          if(_level->onLevel(actlevel)==0) {
594 594
            _el = actlevel;
595 595
            return false;
596 596
          }
597 597
        }
598 598
      next_l:
599 599
        ;
600 600
      }
601 601
      return true;
602 602
    }
603 603

	
604 604
    /// Runs the algorithm.
605 605

	
606 606
    /// This function runs the algorithm.
607 607
    ///
608 608
    /// \return \c true if a feasible circulation is found.
609 609
    ///
610 610
    /// \note Apart from the return value, c.run() is just a shortcut of
611 611
    /// the following code.
612 612
    /// \code
613 613
    ///   c.greedyInit();
614 614
    ///   c.start();
615 615
    /// \endcode
616 616
    bool run() {
617 617
      greedyInit();
618 618
      return start();
619 619
    }
620 620

	
621 621
    /// @}
622 622

	
623 623
    /// \name Query Functions
624 624
    /// The results of the circulation algorithm can be obtained using
625 625
    /// these functions.\n
626 626
    /// Either \ref run() or \ref start() should be called before
627 627
    /// using them.
628 628

	
629 629
    ///@{
630 630

	
631 631
    /// \brief Returns the flow on the given arc.
632 632
    ///
633 633
    /// Returns the flow on the given arc.
634 634
    ///
635 635
    /// \pre Either \ref run() or \ref init() must be called before
636 636
    /// using this function.
637 637
    Value flow(const Arc& arc) const {
638 638
      return (*_flow)[arc];
639 639
    }
640 640

	
641 641
    /// \brief Returns a const reference to the flow map.
642 642
    ///
643 643
    /// Returns a const reference to the arc map storing the found flow.
644 644
    ///
645 645
    /// \pre Either \ref run() or \ref init() must be called before
646 646
    /// using this function.
647
    const FlowMap& flowMap() {
647
    const FlowMap& flowMap() const {
648 648
      return *_flow;
649 649
    }
650 650

	
651 651
    /**
652 652
       \brief Returns \c true if the given node is in a barrier.
653 653

	
654 654
       Barrier is a set \e B of nodes for which
655 655

	
656 656
       \f[ \sum_{a\in\delta_{out}(B)} upper(a) -
657 657
           \sum_{a\in\delta_{in}(B)} lower(a) < \sum_{v\in B}delta(v) \f]
658 658

	
659 659
       holds. The existence of a set with this property prooves that a
660 660
       feasible circualtion cannot exist.
661 661

	
662 662
       This function returns \c true if the given node is in the found
663 663
       barrier. If a feasible circulation is found, the function
664 664
       gives back \c false for every node.
665 665

	
666 666
       \pre Either \ref run() or \ref init() must be called before
667 667
       using this function.
668 668

	
669 669
       \sa barrierMap()
670 670
       \sa checkBarrier()
671 671
    */
672
    bool barrier(const Node& node)
672
    bool barrier(const Node& node) const
673 673
    {
674 674
      return (*_level)[node] >= _el;
675 675
    }
676 676

	
677 677
    /// \brief Gives back a barrier.
678 678
    ///
679 679
    /// This function sets \c bar to the characteristic vector of the
680 680
    /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
681 681
    /// node map with \c bool (or convertible) value type.
682 682
    ///
683 683
    /// If a feasible circulation is found, the function gives back an
684 684
    /// empty set, so \c bar[v] will be \c false for all nodes \c v.
685 685
    ///
686 686
    /// \note This function calls \ref barrier() for each node,
687 687
    /// so it runs in \f$O(n)\f$ time.
688 688
    ///
689 689
    /// \pre Either \ref run() or \ref init() must be called before
690 690
    /// using this function.
691 691
    ///
692 692
    /// \sa barrier()
693 693
    /// \sa checkBarrier()
694 694
    template<class BarrierMap>
695
    void barrierMap(BarrierMap &bar)
695
    void barrierMap(BarrierMap &bar) const
696 696
    {
697 697
      for(NodeIt n(_g);n!=INVALID;++n)
698 698
        bar.set(n, (*_level)[n] >= _el);
699 699
    }
700 700

	
701 701
    /// @}
702 702

	
703 703
    /// \name Checker Functions
704 704
    /// The feasibility of the results can be checked using
705 705
    /// these functions.\n
706 706
    /// Either \ref run() or \ref start() should be called before
707 707
    /// using them.
708 708

	
709 709
    ///@{
710 710

	
711 711
    ///Check if the found flow is a feasible circulation
712 712

	
713 713
    ///Check if the found flow is a feasible circulation,
714 714
    ///
715
    bool checkFlow() {
715
    bool checkFlow() const {
716 716
      for(ArcIt e(_g);e!=INVALID;++e)
717 717
        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
718 718
      for(NodeIt n(_g);n!=INVALID;++n)
719 719
        {
720 720
          Value dif=-(*_delta)[n];
721 721
          for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
722 722
          for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
723 723
          if(_tol.negative(dif)) return false;
724 724
        }
725 725
      return true;
726 726
    }
727 727

	
728 728
    ///Check whether or not the last execution provides a barrier
729 729

	
730 730
    ///Check whether or not the last execution provides a barrier.
731 731
    ///\sa barrier()
732 732
    ///\sa barrierMap()
733
    bool checkBarrier()
733
    bool checkBarrier() const
734 734
    {
735 735
      Value delta=0;
736 736
      for(NodeIt n(_g);n!=INVALID;++n)
737 737
        if(barrier(n))
738 738
          delta-=(*_delta)[n];
739 739
      for(ArcIt e(_g);e!=INVALID;++e)
740 740
        {
741 741
          Node s=_g.source(e);
742 742
          Node t=_g.target(e);
743 743
          if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
744 744
          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
745 745
        }
746 746
      return _tol.negative(delta);
747 747
    }
748 748

	
749 749
    /// @}
750 750

	
751 751
  };
752 752

	
753 753
}
754 754

	
755 755
#endif
Ignore white space 6 line context
... ...
@@ -1116,521 +1116,521 @@
1116 1116
  ///\sa Dfs
1117 1117
  template<class GR>
1118 1118
  DfsWizard<DfsWizardBase<GR> >
1119 1119
  dfs(const GR &digraph)
1120 1120
  {
1121 1121
    return DfsWizard<DfsWizardBase<GR> >(digraph);
1122 1122
  }
1123 1123

	
1124 1124
#ifdef DOXYGEN
1125 1125
  /// \brief Visitor class for DFS.
1126 1126
  ///
1127 1127
  /// This class defines the interface of the DfsVisit events, and
1128 1128
  /// it could be the base of a real visitor class.
1129 1129
  template <typename _Digraph>
1130 1130
  struct DfsVisitor {
1131 1131
    typedef _Digraph Digraph;
1132 1132
    typedef typename Digraph::Arc Arc;
1133 1133
    typedef typename Digraph::Node Node;
1134 1134
    /// \brief Called for the source node of the DFS.
1135 1135
    ///
1136 1136
    /// This function is called for the source node of the DFS.
1137 1137
    void start(const Node& node) {}
1138 1138
    /// \brief Called when the source node is leaved.
1139 1139
    ///
1140 1140
    /// This function is called when the source node is leaved.
1141 1141
    void stop(const Node& node) {}
1142 1142
    /// \brief Called when a node is reached first time.
1143 1143
    ///
1144 1144
    /// This function is called when a node is reached first time.
1145 1145
    void reach(const Node& node) {}
1146 1146
    /// \brief Called when an arc reaches a new node.
1147 1147
    ///
1148 1148
    /// This function is called when the DFS finds an arc whose target node
1149 1149
    /// is not reached yet.
1150 1150
    void discover(const Arc& arc) {}
1151 1151
    /// \brief Called when an arc is examined but its target node is
1152 1152
    /// already discovered.
1153 1153
    ///
1154 1154
    /// This function is called when an arc is examined but its target node is
1155 1155
    /// already discovered.
1156 1156
    void examine(const Arc& arc) {}
1157 1157
    /// \brief Called when the DFS steps back from a node.
1158 1158
    ///
1159 1159
    /// This function is called when the DFS steps back from a node.
1160 1160
    void leave(const Node& node) {}
1161 1161
    /// \brief Called when the DFS steps back on an arc.
1162 1162
    ///
1163 1163
    /// This function is called when the DFS steps back on an arc.
1164 1164
    void backtrack(const Arc& arc) {}
1165 1165
  };
1166 1166
#else
1167 1167
  template <typename _Digraph>
1168 1168
  struct DfsVisitor {
1169 1169
    typedef _Digraph Digraph;
1170 1170
    typedef typename Digraph::Arc Arc;
1171 1171
    typedef typename Digraph::Node Node;
1172 1172
    void start(const Node&) {}
1173 1173
    void stop(const Node&) {}
1174 1174
    void reach(const Node&) {}
1175 1175
    void discover(const Arc&) {}
1176 1176
    void examine(const Arc&) {}
1177 1177
    void leave(const Node&) {}
1178 1178
    void backtrack(const Arc&) {}
1179 1179

	
1180 1180
    template <typename _Visitor>
1181 1181
    struct Constraints {
1182 1182
      void constraints() {
1183 1183
        Arc arc;
1184 1184
        Node node;
1185 1185
        visitor.start(node);
1186 1186
        visitor.stop(arc);
1187 1187
        visitor.reach(node);
1188 1188
        visitor.discover(arc);
1189 1189
        visitor.examine(arc);
1190 1190
        visitor.leave(node);
1191 1191
        visitor.backtrack(arc);
1192 1192
      }
1193 1193
      _Visitor& visitor;
1194 1194
    };
1195 1195
  };
1196 1196
#endif
1197 1197

	
1198 1198
  /// \brief Default traits class of DfsVisit class.
1199 1199
  ///
1200 1200
  /// Default traits class of DfsVisit class.
1201 1201
  /// \tparam _Digraph The type of the digraph the algorithm runs on.
1202 1202
  template<class _Digraph>
1203 1203
  struct DfsVisitDefaultTraits {
1204 1204

	
1205 1205
    /// \brief The type of the digraph the algorithm runs on.
1206 1206
    typedef _Digraph Digraph;
1207 1207

	
1208 1208
    /// \brief The type of the map that indicates which nodes are reached.
1209 1209
    ///
1210 1210
    /// The type of the map that indicates which nodes are reached.
1211 1211
    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
1212 1212
    typedef typename Digraph::template NodeMap<bool> ReachedMap;
1213 1213

	
1214 1214
    /// \brief Instantiates a ReachedMap.
1215 1215
    ///
1216 1216
    /// This function instantiates a ReachedMap.
1217 1217
    /// \param digraph is the digraph, to which
1218 1218
    /// we would like to define the ReachedMap.
1219 1219
    static ReachedMap *createReachedMap(const Digraph &digraph) {
1220 1220
      return new ReachedMap(digraph);
1221 1221
    }
1222 1222

	
1223 1223
  };
1224 1224

	
1225 1225
  /// \ingroup search
1226 1226
  ///
1227 1227
  /// \brief %DFS algorithm class with visitor interface.
1228 1228
  ///
1229 1229
  /// This class provides an efficient implementation of the %DFS algorithm
1230 1230
  /// with visitor interface.
1231 1231
  ///
1232 1232
  /// The %DfsVisit class provides an alternative interface to the Dfs
1233 1233
  /// class. It works with callback mechanism, the DfsVisit object calls
1234 1234
  /// the member functions of the \c Visitor class on every DFS event.
1235 1235
  ///
1236 1236
  /// This interface of the DFS algorithm should be used in special cases
1237 1237
  /// when extra actions have to be performed in connection with certain
1238 1238
  /// events of the DFS algorithm. Otherwise consider to use Dfs or dfs()
1239 1239
  /// instead.
1240 1240
  ///
1241 1241
  /// \tparam _Digraph The type of the digraph the algorithm runs on.
1242 1242
  /// The default value is
1243 1243
  /// \ref ListDigraph. The value of _Digraph is not used directly by
1244 1244
  /// \ref DfsVisit, it is only passed to \ref DfsVisitDefaultTraits.
1245 1245
  /// \tparam _Visitor The Visitor type that is used by the algorithm.
1246 1246
  /// \ref DfsVisitor "DfsVisitor<_Digraph>" is an empty visitor, which
1247 1247
  /// does not observe the DFS events. If you want to observe the DFS
1248 1248
  /// events, you should implement your own visitor class.
1249 1249
  /// \tparam _Traits Traits class to set various data types used by the
1250 1250
  /// algorithm. The default traits class is
1251 1251
  /// \ref DfsVisitDefaultTraits "DfsVisitDefaultTraits<_Digraph>".
1252 1252
  /// See \ref DfsVisitDefaultTraits for the documentation of
1253 1253
  /// a DFS visit traits class.
1254 1254
#ifdef DOXYGEN
1255 1255
  template <typename _Digraph, typename _Visitor, typename _Traits>
1256 1256
#else
1257 1257
  template <typename _Digraph = ListDigraph,
1258 1258
            typename _Visitor = DfsVisitor<_Digraph>,
1259 1259
            typename _Traits = DfsVisitDefaultTraits<_Digraph> >
1260 1260
#endif
1261 1261
  class DfsVisit {
1262 1262
  public:
1263 1263

	
1264 1264
    ///The traits class.
1265 1265
    typedef _Traits Traits;
1266 1266

	
1267 1267
    ///The type of the digraph the algorithm runs on.
1268 1268
    typedef typename Traits::Digraph Digraph;
1269 1269

	
1270 1270
    ///The visitor type used by the algorithm.
1271 1271
    typedef _Visitor Visitor;
1272 1272

	
1273 1273
    ///The type of the map that indicates which nodes are reached.
1274 1274
    typedef typename Traits::ReachedMap ReachedMap;
1275 1275

	
1276 1276
  private:
1277 1277

	
1278 1278
    typedef typename Digraph::Node Node;
1279 1279
    typedef typename Digraph::NodeIt NodeIt;
1280 1280
    typedef typename Digraph::Arc Arc;
1281 1281
    typedef typename Digraph::OutArcIt OutArcIt;
1282 1282

	
1283 1283
    //Pointer to the underlying digraph.
1284 1284
    const Digraph *_digraph;
1285 1285
    //Pointer to the visitor object.
1286 1286
    Visitor *_visitor;
1287 1287
    //Pointer to the map of reached status of the nodes.
1288 1288
    ReachedMap *_reached;
1289 1289
    //Indicates if _reached is locally allocated (true) or not.
1290 1290
    bool local_reached;
1291 1291

	
1292 1292
    std::vector<typename Digraph::Arc> _stack;
1293 1293
    int _stack_head;
1294 1294

	
1295 1295
    //Creates the maps if necessary.
1296 1296
    void create_maps() {
1297 1297
      if(!_reached) {
1298 1298
        local_reached = true;
1299 1299
        _reached = Traits::createReachedMap(*_digraph);
1300 1300
      }
1301 1301
    }
1302 1302

	
1303 1303
  protected:
1304 1304

	
1305 1305
    DfsVisit() {}
1306 1306

	
1307 1307
  public:
1308 1308

	
1309 1309
    typedef DfsVisit Create;
1310 1310

	
1311 1311
    /// \name Named Template Parameters
1312 1312

	
1313 1313
    ///@{
1314 1314
    template <class T>
1315 1315
    struct SetReachedMapTraits : public Traits {
1316 1316
      typedef T ReachedMap;
1317 1317
      static ReachedMap *createReachedMap(const Digraph &digraph) {
1318 1318
        LEMON_ASSERT(false, "ReachedMap is not initialized");
1319 1319
        return 0; // ignore warnings
1320 1320
      }
1321 1321
    };
1322 1322
    /// \brief \ref named-templ-param "Named parameter" for setting
1323 1323
    /// ReachedMap type.
1324 1324
    ///
1325 1325
    /// \ref named-templ-param "Named parameter" for setting ReachedMap type.
1326 1326
    template <class T>
1327 1327
    struct SetReachedMap : public DfsVisit< Digraph, Visitor,
1328 1328
                                            SetReachedMapTraits<T> > {
1329 1329
      typedef DfsVisit< Digraph, Visitor, SetReachedMapTraits<T> > Create;
1330 1330
    };
1331 1331
    ///@}
1332 1332

	
1333 1333
  public:
1334 1334

	
1335 1335
    /// \brief Constructor.
1336 1336
    ///
1337 1337
    /// Constructor.
1338 1338
    ///
1339 1339
    /// \param digraph The digraph the algorithm runs on.
1340 1340
    /// \param visitor The visitor object of the algorithm.
1341 1341
    DfsVisit(const Digraph& digraph, Visitor& visitor)
1342 1342
      : _digraph(&digraph), _visitor(&visitor),
1343 1343
        _reached(0), local_reached(false) {}
1344 1344

	
1345 1345
    /// \brief Destructor.
1346 1346
    ~DfsVisit() {
1347 1347
      if(local_reached) delete _reached;
1348 1348
    }
1349 1349

	
1350 1350
    /// \brief Sets the map that indicates which nodes are reached.
1351 1351
    ///
1352 1352
    /// Sets the map that indicates which nodes are reached.
1353 1353
    /// If you don't use this function before calling \ref run(Node) "run()"
1354 1354
    /// or \ref init(), an instance will be allocated automatically.
1355 1355
    /// The destructor deallocates this automatically allocated map,
1356 1356
    /// of course.
1357 1357
    /// \return <tt> (*this) </tt>
1358 1358
    DfsVisit &reachedMap(ReachedMap &m) {
1359 1359
      if(local_reached) {
1360 1360
        delete _reached;
1361 1361
        local_reached=false;
1362 1362
      }
1363 1363
      _reached = &m;
1364 1364
      return *this;
1365 1365
    }
1366 1366

	
1367 1367
  public:
1368 1368

	
1369 1369
    /// \name Execution Control
1370 1370
    /// The simplest way to execute the DFS algorithm is to use one of the
1371 1371
    /// member functions called \ref run(Node) "run()".\n
1372 1372
    /// If you need more control on the execution, first you have to call
1373 1373
    /// \ref init(), then you can add a source node with \ref addSource()
1374 1374
    /// and perform the actual computation with \ref start().
1375 1375
    /// This procedure can be repeated if there are nodes that have not
1376 1376
    /// been reached.
1377 1377

	
1378 1378
    /// @{
1379 1379

	
1380 1380
    /// \brief Initializes the internal data structures.
1381 1381
    ///
1382 1382
    /// Initializes the internal data structures.
1383 1383
    void init() {
1384 1384
      create_maps();
1385 1385
      _stack.resize(countNodes(*_digraph));
1386 1386
      _stack_head = -1;
1387 1387
      for (NodeIt u(*_digraph) ; u != INVALID ; ++u) {
1388 1388
        _reached->set(u, false);
1389 1389
      }
1390 1390
    }
1391 1391

	
1392 1392
    /// \brief Adds a new source node.
1393 1393
    ///
1394 1394
    /// Adds a new source node to the set of nodes to be processed.
1395 1395
    ///
1396 1396
    /// \pre The stack must be empty. Otherwise the algorithm gives
1397 1397
    /// wrong results. (One of the outgoing arcs of all the source nodes
1398 1398
    /// except for the last one will not be visited and distances will
1399 1399
    /// also be wrong.)
1400 1400
    void addSource(Node s)
1401 1401
    {
1402 1402
      LEMON_DEBUG(emptyQueue(), "The stack is not empty.");
1403 1403
      if(!(*_reached)[s]) {
1404 1404
          _reached->set(s,true);
1405 1405
          _visitor->start(s);
1406 1406
          _visitor->reach(s);
1407 1407
          Arc e;
1408 1408
          _digraph->firstOut(e, s);
1409 1409
          if (e != INVALID) {
1410 1410
            _stack[++_stack_head] = e;
1411 1411
          } else {
1412 1412
            _visitor->leave(s);
1413 1413
          }
1414 1414
        }
1415 1415
    }
1416 1416

	
1417 1417
    /// \brief Processes the next arc.
1418 1418
    ///
1419 1419
    /// Processes the next arc.
1420 1420
    ///
1421 1421
    /// \return The processed arc.
1422 1422
    ///
1423 1423
    /// \pre The stack must not be empty.
1424 1424
    Arc processNextArc() {
1425 1425
      Arc e = _stack[_stack_head];
1426 1426
      Node m = _digraph->target(e);
1427 1427
      if(!(*_reached)[m]) {
1428 1428
        _visitor->discover(e);
1429 1429
        _visitor->reach(m);
1430 1430
        _reached->set(m, true);
1431 1431
        _digraph->firstOut(_stack[++_stack_head], m);
1432 1432
      } else {
1433 1433
        _visitor->examine(e);
1434 1434
        m = _digraph->source(e);
1435 1435
        _digraph->nextOut(_stack[_stack_head]);
1436 1436
      }
1437 1437
      while (_stack_head>=0 && _stack[_stack_head] == INVALID) {
1438 1438
        _visitor->leave(m);
1439 1439
        --_stack_head;
1440 1440
        if (_stack_head >= 0) {
1441 1441
          _visitor->backtrack(_stack[_stack_head]);
1442 1442
          m = _digraph->source(_stack[_stack_head]);
1443 1443
          _digraph->nextOut(_stack[_stack_head]);
1444 1444
        } else {
1445 1445
          _visitor->stop(m);
1446 1446
        }
1447 1447
      }
1448 1448
      return e;
1449 1449
    }
1450 1450

	
1451 1451
    /// \brief Next arc to be processed.
1452 1452
    ///
1453 1453
    /// Next arc to be processed.
1454 1454
    ///
1455 1455
    /// \return The next arc to be processed or INVALID if the stack is
1456 1456
    /// empty.
1457 1457
    Arc nextArc() const {
1458 1458
      return _stack_head >= 0 ? _stack[_stack_head] : INVALID;
1459 1459
    }
1460 1460

	
1461 1461
    /// \brief Returns \c false if there are nodes
1462 1462
    /// to be processed.
1463 1463
    ///
1464 1464
    /// Returns \c false if there are nodes
1465 1465
    /// to be processed in the queue (stack).
1466 1466
    bool emptyQueue() const { return _stack_head < 0; }
1467 1467

	
1468 1468
    /// \brief Returns the number of the nodes to be processed.
1469 1469
    ///
1470 1470
    /// Returns the number of the nodes to be processed in the queue (stack).
1471 1471
    int queueSize() const { return _stack_head + 1; }
1472 1472

	
1473 1473
    /// \brief Executes the algorithm.
1474 1474
    ///
1475 1475
    /// Executes the algorithm.
1476 1476
    ///
1477 1477
    /// This method runs the %DFS algorithm from the root node
1478 1478
    /// in order to compute the %DFS path to each node.
1479 1479
    ///
1480 1480
    /// The algorithm computes
1481 1481
    /// - the %DFS tree,
1482 1482
    /// - the distance of each node from the root in the %DFS tree.
1483 1483
    ///
1484 1484
    /// \pre init() must be called and a root node should be
1485 1485
    /// added with addSource() before using this function.
1486 1486
    ///
1487 1487
    /// \note <tt>d.start()</tt> is just a shortcut of the following code.
1488 1488
    /// \code
1489 1489
    ///   while ( !d.emptyQueue() ) {
1490 1490
    ///     d.processNextArc();
1491 1491
    ///   }
1492 1492
    /// \endcode
1493 1493
    void start() {
1494 1494
      while ( !emptyQueue() ) processNextArc();
1495 1495
    }
1496 1496

	
1497 1497
    /// \brief Executes the algorithm until the given target node is reached.
1498 1498
    ///
1499 1499
    /// Executes the algorithm until the given target node is reached.
1500 1500
    ///
1501 1501
    /// This method runs the %DFS algorithm from the root node
1502 1502
    /// in order to compute the DFS path to \c t.
1503 1503
    ///
1504 1504
    /// The algorithm computes
1505 1505
    /// - the %DFS path to \c t,
1506 1506
    /// - the distance of \c t from the root in the %DFS tree.
1507 1507
    ///
1508 1508
    /// \pre init() must be called and a root node should be added
1509 1509
    /// with addSource() before using this function.
1510 1510
    void start(Node t) {
1511 1511
      while ( !emptyQueue() && _digraph->target(_stack[_stack_head]) != t )
1512 1512
        processNextArc();
1513 1513
    }
1514 1514

	
1515 1515
    /// \brief Executes the algorithm until a condition is met.
1516 1516
    ///
1517 1517
    /// Executes the algorithm until a condition is met.
1518 1518
    ///
1519 1519
    /// This method runs the %DFS algorithm from the root node
1520 1520
    /// until an arc \c a with <tt>am[a]</tt> true is found.
1521 1521
    ///
1522 1522
    /// \param am A \c bool (or convertible) arc map. The algorithm
1523 1523
    /// will stop when it reaches an arc \c a with <tt>am[a]</tt> true.
1524 1524
    ///
1525 1525
    /// \return The reached arc \c a with <tt>am[a]</tt> true or
1526 1526
    /// \c INVALID if no such arc was found.
1527 1527
    ///
1528 1528
    /// \pre init() must be called and a root node should be added
1529 1529
    /// with addSource() before using this function.
1530 1530
    ///
1531 1531
    /// \warning Contrary to \ref Bfs and \ref Dijkstra, \c am is an arc map,
1532 1532
    /// not a node map.
1533 1533
    template <typename AM>
1534 1534
    Arc start(const AM &am) {
1535 1535
      while ( !emptyQueue() && !am[_stack[_stack_head]] )
1536 1536
        processNextArc();
1537 1537
      return emptyQueue() ? INVALID : _stack[_stack_head];
1538 1538
    }
1539 1539

	
1540 1540
    /// \brief Runs the algorithm from the given source node.
1541 1541
    ///
1542 1542
    /// This method runs the %DFS algorithm from node \c s.
1543 1543
    /// in order to compute the DFS path to each node.
1544 1544
    ///
1545 1545
    /// The algorithm computes
1546 1546
    /// - the %DFS tree,
1547 1547
    /// - the distance of each node from the root in the %DFS tree.
1548 1548
    ///
1549 1549
    /// \note <tt>d.run(s)</tt> is just a shortcut of the following code.
1550 1550
    ///\code
1551 1551
    ///   d.init();
1552 1552
    ///   d.addSource(s);
1553 1553
    ///   d.start();
1554 1554
    ///\endcode
1555 1555
    void run(Node s) {
1556 1556
      init();
1557 1557
      addSource(s);
1558 1558
      start();
1559 1559
    }
1560 1560

	
1561 1561
    /// \brief Finds the %DFS path between \c s and \c t.
1562 1562

	
1563 1563
    /// This method runs the %DFS algorithm from node \c s
1564 1564
    /// in order to compute the DFS path to node \c t
1565 1565
    /// (it stops searching when \c t is processed).
1566 1566
    ///
1567 1567
    /// \return \c true if \c t is reachable form \c s.
1568 1568
    ///
1569 1569
    /// \note Apart from the return value, <tt>d.run(s,t)</tt> is
1570 1570
    /// just a shortcut of the following code.
1571 1571
    ///\code
1572 1572
    ///   d.init();
1573 1573
    ///   d.addSource(s);
1574 1574
    ///   d.start(t);
1575 1575
    ///\endcode
1576 1576
    bool run(Node s,Node t) {
1577 1577
      init();
1578 1578
      addSource(s);
1579 1579
      start(t);
1580 1580
      return reached(t);
1581 1581
    }
1582 1582

	
1583 1583
    /// \brief Runs the algorithm to visit all nodes in the digraph.
1584 1584

	
1585 1585
    /// This method runs the %DFS algorithm in order to
1586 1586
    /// compute the %DFS path to each node.
1587 1587
    ///
1588 1588
    /// The algorithm computes
1589 1589
    /// - the %DFS tree (forest),
1590 1590
    /// - the distance of each node from the root(s) in the %DFS tree.
1591 1591
    ///
1592 1592
    /// \note <tt>d.run()</tt> is just a shortcut of the following code.
1593 1593
    ///\code
1594 1594
    ///   d.init();
1595 1595
    ///   for (NodeIt n(digraph); n != INVALID; ++n) {
1596 1596
    ///     if (!d.reached(n)) {
1597 1597
    ///       d.addSource(n);
1598 1598
    ///       d.start();
1599 1599
    ///     }
1600 1600
    ///   }
1601 1601
    ///\endcode
1602 1602
    void run() {
1603 1603
      init();
1604 1604
      for (NodeIt it(*_digraph); it != INVALID; ++it) {
1605 1605
        if (!reached(it)) {
1606 1606
          addSource(it);
1607 1607
          start();
1608 1608
        }
1609 1609
      }
1610 1610
    }
1611 1611

	
1612 1612
    ///@}
1613 1613

	
1614 1614
    /// \name Query Functions
1615 1615
    /// The results of the DFS algorithm can be obtained using these
1616 1616
    /// functions.\n
1617 1617
    /// Either \ref run(Node) "run()" or \ref start() should be called
1618 1618
    /// before using them.
1619 1619

	
1620 1620
    ///@{
1621 1621

	
1622 1622
    /// \brief Checks if a node is reached from the root(s).
1623 1623
    ///
1624 1624
    /// Returns \c true if \c v is reached from the root(s).
1625 1625
    ///
1626 1626
    /// \pre Either \ref run(Node) "run()" or \ref init()
1627 1627
    /// must be called before using this function.
1628
    bool reached(Node v) { return (*_reached)[v]; }
1628
    bool reached(Node v) const { return (*_reached)[v]; }
1629 1629

	
1630 1630
    ///@}
1631 1631

	
1632 1632
  };
1633 1633

	
1634 1634
} //END OF NAMESPACE LEMON
1635 1635

	
1636 1636
#endif
Ignore white space 6 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_PREFLOW_H
20 20
#define LEMON_PREFLOW_H
21 21

	
22 22
#include <lemon/tolerance.h>
23 23
#include <lemon/elevator.h>
24 24

	
25 25
/// \file
26 26
/// \ingroup max_flow
27 27
/// \brief Implementation of the preflow algorithm.
28 28

	
29 29
namespace lemon {
30 30

	
31 31
  /// \brief Default traits class of Preflow class.
32 32
  ///
33 33
  /// Default traits class of Preflow class.
34 34
  /// \tparam _Digraph Digraph type.
35 35
  /// \tparam _CapacityMap Capacity map type.
36 36
  template <typename _Digraph, typename _CapacityMap>
37 37
  struct PreflowDefaultTraits {
38 38

	
39 39
    /// \brief The type of the digraph the algorithm runs on.
40 40
    typedef _Digraph Digraph;
41 41

	
42 42
    /// \brief The type of the map that stores the arc capacities.
43 43
    ///
44 44
    /// The type of the map that stores the arc capacities.
45 45
    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
46 46
    typedef _CapacityMap CapacityMap;
47 47

	
48 48
    /// \brief The type of the flow values.
49 49
    typedef typename CapacityMap::Value Value;
50 50

	
51 51
    /// \brief The type of the map that stores the flow values.
52 52
    ///
53 53
    /// The type of the map that stores the flow values.
54 54
    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
55 55
    typedef typename Digraph::template ArcMap<Value> FlowMap;
56 56

	
57 57
    /// \brief Instantiates a FlowMap.
58 58
    ///
59 59
    /// This function instantiates a \ref FlowMap.
60 60
    /// \param digraph The digraph, to which we would like to define
61 61
    /// the flow map.
62 62
    static FlowMap* createFlowMap(const Digraph& digraph) {
63 63
      return new FlowMap(digraph);
64 64
    }
65 65

	
66 66
    /// \brief The elevator type used by Preflow algorithm.
67 67
    ///
68 68
    /// The elevator type used by Preflow algorithm.
69 69
    ///
70 70
    /// \sa Elevator
71 71
    /// \sa LinkedElevator
72 72
    typedef LinkedElevator<Digraph, typename Digraph::Node> Elevator;
73 73

	
74 74
    /// \brief Instantiates an Elevator.
75 75
    ///
76 76
    /// This function instantiates an \ref Elevator.
77 77
    /// \param digraph The digraph, to which we would like to define
78 78
    /// the elevator.
79 79
    /// \param max_level The maximum level of the elevator.
80 80
    static Elevator* createElevator(const Digraph& digraph, int max_level) {
81 81
      return new Elevator(digraph, max_level);
82 82
    }
83 83

	
84 84
    /// \brief The tolerance used by the algorithm
85 85
    ///
86 86
    /// The tolerance used by the algorithm to handle inexact computation.
87 87
    typedef lemon::Tolerance<Value> Tolerance;
88 88

	
89 89
  };
90 90

	
91 91

	
92 92
  /// \ingroup max_flow
93 93
  ///
94 94
  /// \brief %Preflow algorithm class.
95 95
  ///
96 96
  /// This class provides an implementation of Goldberg-Tarjan's \e preflow
97 97
  /// \e push-relabel algorithm producing a flow of maximum value in a
98 98
  /// digraph. The preflow algorithms are the fastest known maximum
99 99
  /// flow algorithms. The current implementation use a mixture of the
100 100
  /// \e "highest label" and the \e "bound decrease" heuristics.
101 101
  /// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$.
102 102
  ///
103 103
  /// The algorithm consists of two phases. After the first phase
104 104
  /// the maximum flow value and the minimum cut is obtained. The
105 105
  /// second phase constructs a feasible maximum flow on each arc.
106 106
  ///
107 107
  /// \tparam _Digraph The type of the digraph the algorithm runs on.
108 108
  /// \tparam _CapacityMap The type of the capacity map. The default map
109 109
  /// type is \ref concepts::Digraph::ArcMap "_Digraph::ArcMap<int>".
110 110
#ifdef DOXYGEN
111 111
  template <typename _Digraph, typename _CapacityMap, typename _Traits>
112 112
#else
113 113
  template <typename _Digraph,
114 114
            typename _CapacityMap = typename _Digraph::template ArcMap<int>,
115 115
            typename _Traits = PreflowDefaultTraits<_Digraph, _CapacityMap> >
116 116
#endif
117 117
  class Preflow {
118 118
  public:
119 119

	
120 120
    ///The \ref PreflowDefaultTraits "traits class" of the algorithm.
121 121
    typedef _Traits Traits;
122 122
    ///The type of the digraph the algorithm runs on.
123 123
    typedef typename Traits::Digraph Digraph;
124 124
    ///The type of the capacity map.
125 125
    typedef typename Traits::CapacityMap CapacityMap;
126 126
    ///The type of the flow values.
127 127
    typedef typename Traits::Value Value;
128 128

	
129 129
    ///The type of the flow map.
130 130
    typedef typename Traits::FlowMap FlowMap;
131 131
    ///The type of the elevator.
132 132
    typedef typename Traits::Elevator Elevator;
133 133
    ///The type of the tolerance.
134 134
    typedef typename Traits::Tolerance Tolerance;
135 135

	
136 136
  private:
137 137

	
138 138
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
139 139

	
140 140
    const Digraph& _graph;
141 141
    const CapacityMap* _capacity;
142 142

	
143 143
    int _node_num;
144 144

	
145 145
    Node _source, _target;
146 146

	
147 147
    FlowMap* _flow;
148 148
    bool _local_flow;
149 149

	
150 150
    Elevator* _level;
151 151
    bool _local_level;
152 152

	
153 153
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
154 154
    ExcessMap* _excess;
155 155

	
156 156
    Tolerance _tolerance;
157 157

	
158 158
    bool _phase;
159 159

	
160 160

	
161 161
    void createStructures() {
162 162
      _node_num = countNodes(_graph);
163 163

	
164 164
      if (!_flow) {
165 165
        _flow = Traits::createFlowMap(_graph);
166 166
        _local_flow = true;
167 167
      }
168 168
      if (!_level) {
169 169
        _level = Traits::createElevator(_graph, _node_num);
170 170
        _local_level = true;
171 171
      }
172 172
      if (!_excess) {
173 173
        _excess = new ExcessMap(_graph);
174 174
      }
175 175
    }
176 176

	
177 177
    void destroyStructures() {
178 178
      if (_local_flow) {
179 179
        delete _flow;
180 180
      }
181 181
      if (_local_level) {
182 182
        delete _level;
183 183
      }
184 184
      if (_excess) {
185 185
        delete _excess;
186 186
      }
187 187
    }
188 188

	
189 189
  public:
190 190

	
191 191
    typedef Preflow Create;
192 192

	
193 193
    ///\name Named Template Parameters
194 194

	
195 195
    ///@{
196 196

	
197 197
    template <typename _FlowMap>
198 198
    struct SetFlowMapTraits : public Traits {
199 199
      typedef _FlowMap FlowMap;
200 200
      static FlowMap *createFlowMap(const Digraph&) {
201 201
        LEMON_ASSERT(false, "FlowMap is not initialized");
202 202
        return 0; // ignore warnings
203 203
      }
204 204
    };
205 205

	
206 206
    /// \brief \ref named-templ-param "Named parameter" for setting
207 207
    /// FlowMap type
208 208
    ///
209 209
    /// \ref named-templ-param "Named parameter" for setting FlowMap
210 210
    /// type.
211 211
    template <typename _FlowMap>
212 212
    struct SetFlowMap
213 213
      : public Preflow<Digraph, CapacityMap, SetFlowMapTraits<_FlowMap> > {
214 214
      typedef Preflow<Digraph, CapacityMap,
215 215
                      SetFlowMapTraits<_FlowMap> > Create;
216 216
    };
217 217

	
218 218
    template <typename _Elevator>
219 219
    struct SetElevatorTraits : public Traits {
220 220
      typedef _Elevator Elevator;
221 221
      static Elevator *createElevator(const Digraph&, int) {
222 222
        LEMON_ASSERT(false, "Elevator is not initialized");
223 223
        return 0; // ignore warnings
224 224
      }
225 225
    };
226 226

	
227 227
    /// \brief \ref named-templ-param "Named parameter" for setting
228 228
    /// Elevator type
229 229
    ///
230 230
    /// \ref named-templ-param "Named parameter" for setting Elevator
231 231
    /// type. If this named parameter is used, then an external
232 232
    /// elevator object must be passed to the algorithm using the
233 233
    /// \ref elevator(Elevator&) "elevator()" function before calling
234 234
    /// \ref run() or \ref init().
235 235
    /// \sa SetStandardElevator
236 236
    template <typename _Elevator>
237 237
    struct SetElevator
238 238
      : public Preflow<Digraph, CapacityMap, SetElevatorTraits<_Elevator> > {
239 239
      typedef Preflow<Digraph, CapacityMap,
240 240
                      SetElevatorTraits<_Elevator> > Create;
241 241
    };
242 242

	
243 243
    template <typename _Elevator>
244 244
    struct SetStandardElevatorTraits : public Traits {
245 245
      typedef _Elevator Elevator;
246 246
      static Elevator *createElevator(const Digraph& digraph, int max_level) {
247 247
        return new Elevator(digraph, max_level);
248 248
      }
249 249
    };
250 250

	
251 251
    /// \brief \ref named-templ-param "Named parameter" for setting
252 252
    /// Elevator type with automatic allocation
253 253
    ///
254 254
    /// \ref named-templ-param "Named parameter" for setting Elevator
255 255
    /// type with automatic allocation.
256 256
    /// The Elevator should have standard constructor interface to be
257 257
    /// able to automatically created by the algorithm (i.e. the
258 258
    /// digraph and the maximum level should be passed to it).
259 259
    /// However an external elevator object could also be passed to the
260 260
    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
261 261
    /// before calling \ref run() or \ref init().
262 262
    /// \sa SetElevator
263 263
    template <typename _Elevator>
264 264
    struct SetStandardElevator
265 265
      : public Preflow<Digraph, CapacityMap,
266 266
                       SetStandardElevatorTraits<_Elevator> > {
267 267
      typedef Preflow<Digraph, CapacityMap,
268 268
                      SetStandardElevatorTraits<_Elevator> > Create;
269 269
    };
270 270

	
271 271
    /// @}
272 272

	
273 273
  protected:
274 274

	
275 275
    Preflow() {}
276 276

	
277 277
  public:
278 278

	
279 279

	
280 280
    /// \brief The constructor of the class.
281 281
    ///
282 282
    /// The constructor of the class.
283 283
    /// \param digraph The digraph the algorithm runs on.
284 284
    /// \param capacity The capacity of the arcs.
285 285
    /// \param source The source node.
286 286
    /// \param target The target node.
287 287
    Preflow(const Digraph& digraph, const CapacityMap& capacity,
288 288
            Node source, Node target)
289 289
      : _graph(digraph), _capacity(&capacity),
290 290
        _node_num(0), _source(source), _target(target),
291 291
        _flow(0), _local_flow(false),
292 292
        _level(0), _local_level(false),
293 293
        _excess(0), _tolerance(), _phase() {}
294 294

	
295 295
    /// \brief Destructor.
296 296
    ///
297 297
    /// Destructor.
298 298
    ~Preflow() {
299 299
      destroyStructures();
300 300
    }
301 301

	
302 302
    /// \brief Sets the capacity map.
303 303
    ///
304 304
    /// Sets the capacity map.
305 305
    /// \return <tt>(*this)</tt>
306 306
    Preflow& capacityMap(const CapacityMap& map) {
307 307
      _capacity = &map;
308 308
      return *this;
309 309
    }
310 310

	
311 311
    /// \brief Sets the flow map.
312 312
    ///
313 313
    /// Sets the flow map.
314 314
    /// If you don't use this function before calling \ref run() or
315 315
    /// \ref init(), an instance will be allocated automatically.
316 316
    /// The destructor deallocates this automatically allocated map,
317 317
    /// of course.
318 318
    /// \return <tt>(*this)</tt>
319 319
    Preflow& flowMap(FlowMap& map) {
320 320
      if (_local_flow) {
321 321
        delete _flow;
322 322
        _local_flow = false;
323 323
      }
324 324
      _flow = &map;
325 325
      return *this;
326 326
    }
327 327

	
328 328
    /// \brief Sets the source node.
329 329
    ///
330 330
    /// Sets the source node.
331 331
    /// \return <tt>(*this)</tt>
332 332
    Preflow& source(const Node& node) {
333 333
      _source = node;
334 334
      return *this;
335 335
    }
336 336

	
337 337
    /// \brief Sets the target node.
338 338
    ///
339 339
    /// Sets the target node.
340 340
    /// \return <tt>(*this)</tt>
341 341
    Preflow& target(const Node& node) {
342 342
      _target = node;
343 343
      return *this;
344 344
    }
345 345

	
346 346
    /// \brief Sets the elevator used by algorithm.
347 347
    ///
348 348
    /// Sets the elevator used by algorithm.
349 349
    /// If you don't use this function before calling \ref run() or
350 350
    /// \ref init(), an instance will be allocated automatically.
351 351
    /// The destructor deallocates this automatically allocated elevator,
352 352
    /// of course.
353 353
    /// \return <tt>(*this)</tt>
354 354
    Preflow& elevator(Elevator& elevator) {
355 355
      if (_local_level) {
356 356
        delete _level;
357 357
        _local_level = false;
358 358
      }
359 359
      _level = &elevator;
360 360
      return *this;
361 361
    }
362 362

	
363 363
    /// \brief Returns a const reference to the elevator.
364 364
    ///
365 365
    /// Returns a const reference to the elevator.
366 366
    ///
367 367
    /// \pre Either \ref run() or \ref init() must be called before
368 368
    /// using this function.
369
    const Elevator& elevator() {
369
    const Elevator& elevator() const {
370 370
      return *_level;
371 371
    }
372 372

	
373 373
    /// \brief Sets the tolerance used by algorithm.
374 374
    ///
375 375
    /// Sets the tolerance used by algorithm.
376 376
    Preflow& tolerance(const Tolerance& tolerance) const {
377 377
      _tolerance = tolerance;
378 378
      return *this;
379 379
    }
380 380

	
381 381
    /// \brief Returns a const reference to the tolerance.
382 382
    ///
383 383
    /// Returns a const reference to the tolerance.
384 384
    const Tolerance& tolerance() const {
385 385
      return tolerance;
386 386
    }
387 387

	
388 388
    /// \name Execution Control
389 389
    /// The simplest way to execute the preflow algorithm is to use
390 390
    /// \ref run() or \ref runMinCut().\n
391 391
    /// If you need more control on the initial solution or the execution,
392 392
    /// first you have to call one of the \ref init() functions, then
393 393
    /// \ref startFirstPhase() and if you need it \ref startSecondPhase().
394 394

	
395 395
    ///@{
396 396

	
397 397
    /// \brief Initializes the internal data structures.
398 398
    ///
399 399
    /// Initializes the internal data structures and sets the initial
400 400
    /// flow to zero on each arc.
401 401
    void init() {
402 402
      createStructures();
403 403

	
404 404
      _phase = true;
405 405
      for (NodeIt n(_graph); n != INVALID; ++n) {
406 406
        _excess->set(n, 0);
407 407
      }
408 408

	
409 409
      for (ArcIt e(_graph); e != INVALID; ++e) {
410 410
        _flow->set(e, 0);
411 411
      }
412 412

	
413 413
      typename Digraph::template NodeMap<bool> reached(_graph, false);
414 414

	
415 415
      _level->initStart();
416 416
      _level->initAddItem(_target);
417 417

	
418 418
      std::vector<Node> queue;
419 419
      reached.set(_source, true);
420 420

	
421 421
      queue.push_back(_target);
422 422
      reached.set(_target, true);
423 423
      while (!queue.empty()) {
424 424
        _level->initNewLevel();
425 425
        std::vector<Node> nqueue;
426 426
        for (int i = 0; i < int(queue.size()); ++i) {
427 427
          Node n = queue[i];
428 428
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
429 429
            Node u = _graph.source(e);
430 430
            if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
431 431
              reached.set(u, true);
432 432
              _level->initAddItem(u);
433 433
              nqueue.push_back(u);
434 434
            }
435 435
          }
436 436
        }
437 437
        queue.swap(nqueue);
438 438
      }
439 439
      _level->initFinish();
440 440

	
441 441
      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
442 442
        if (_tolerance.positive((*_capacity)[e])) {
443 443
          Node u = _graph.target(e);
444 444
          if ((*_level)[u] == _level->maxLevel()) continue;
445 445
          _flow->set(e, (*_capacity)[e]);
446 446
          _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
447 447
          if (u != _target && !_level->active(u)) {
448 448
            _level->activate(u);
449 449
          }
450 450
        }
451 451
      }
452 452
    }
453 453

	
454 454
    /// \brief Initializes the internal data structures using the
455 455
    /// given flow map.
456 456
    ///
457 457
    /// Initializes the internal data structures and sets the initial
458 458
    /// flow to the given \c flowMap. The \c flowMap should contain a
459 459
    /// flow or at least a preflow, i.e. at each node excluding the
460 460
    /// source node the incoming flow should greater or equal to the
461 461
    /// outgoing flow.
462 462
    /// \return \c false if the given \c flowMap is not a preflow.
463 463
    template <typename FlowMap>
464 464
    bool init(const FlowMap& flowMap) {
465 465
      createStructures();
466 466

	
467 467
      for (ArcIt e(_graph); e != INVALID; ++e) {
468 468
        _flow->set(e, flowMap[e]);
469 469
      }
470 470

	
471 471
      for (NodeIt n(_graph); n != INVALID; ++n) {
472 472
        Value excess = 0;
473 473
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
474 474
          excess += (*_flow)[e];
475 475
        }
476 476
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
477 477
          excess -= (*_flow)[e];
478 478
        }
479 479
        if (excess < 0 && n != _source) return false;
480 480
        _excess->set(n, excess);
481 481
      }
482 482

	
483 483
      typename Digraph::template NodeMap<bool> reached(_graph, false);
484 484

	
485 485
      _level->initStart();
486 486
      _level->initAddItem(_target);
487 487

	
488 488
      std::vector<Node> queue;
489 489
      reached.set(_source, true);
490 490

	
491 491
      queue.push_back(_target);
492 492
      reached.set(_target, true);
493 493
      while (!queue.empty()) {
494 494
        _level->initNewLevel();
495 495
        std::vector<Node> nqueue;
496 496
        for (int i = 0; i < int(queue.size()); ++i) {
497 497
          Node n = queue[i];
498 498
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
499 499
            Node u = _graph.source(e);
500 500
            if (!reached[u] &&
501 501
                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
502 502
              reached.set(u, true);
503 503
              _level->initAddItem(u);
504 504
              nqueue.push_back(u);
505 505
            }
506 506
          }
507 507
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
508 508
            Node v = _graph.target(e);
509 509
            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
510 510
              reached.set(v, true);
511 511
              _level->initAddItem(v);
512 512
              nqueue.push_back(v);
513 513
            }
514 514
          }
515 515
        }
516 516
        queue.swap(nqueue);
517 517
      }
518 518
      _level->initFinish();
519 519

	
520 520
      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
521 521
        Value rem = (*_capacity)[e] - (*_flow)[e];
522 522
        if (_tolerance.positive(rem)) {
523 523
          Node u = _graph.target(e);
524 524
          if ((*_level)[u] == _level->maxLevel()) continue;
525 525
          _flow->set(e, (*_capacity)[e]);
526 526
          _excess->set(u, (*_excess)[u] + rem);
527 527
          if (u != _target && !_level->active(u)) {
528 528
            _level->activate(u);
529 529
          }
530 530
        }
531 531
      }
532 532
      for (InArcIt e(_graph, _source); e != INVALID; ++e) {
533 533
        Value rem = (*_flow)[e];
534 534
        if (_tolerance.positive(rem)) {
535 535
          Node v = _graph.source(e);
536 536
          if ((*_level)[v] == _level->maxLevel()) continue;
537 537
          _flow->set(e, 0);
538 538
          _excess->set(v, (*_excess)[v] + rem);
539 539
          if (v != _target && !_level->active(v)) {
540 540
            _level->activate(v);
541 541
          }
542 542
        }
543 543
      }
544 544
      return true;
545 545
    }
546 546

	
547 547
    /// \brief Starts the first phase of the preflow algorithm.
548 548
    ///
549 549
    /// The preflow algorithm consists of two phases, this method runs
550 550
    /// the first phase. After the first phase the maximum flow value
551 551
    /// and a minimum value cut can already be computed, although a
552 552
    /// maximum flow is not yet obtained. So after calling this method
553 553
    /// \ref flowValue() returns the value of a maximum flow and \ref
554 554
    /// minCut() returns a minimum cut.
555 555
    /// \pre One of the \ref init() functions must be called before
556 556
    /// using this function.
557 557
    void startFirstPhase() {
558 558
      _phase = true;
559 559

	
560 560
      Node n = _level->highestActive();
561 561
      int level = _level->highestActiveLevel();
562 562
      while (n != INVALID) {
563 563
        int num = _node_num;
564 564

	
565 565
        while (num > 0 && n != INVALID) {
566 566
          Value excess = (*_excess)[n];
567 567
          int new_level = _level->maxLevel();
568 568

	
569 569
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
570 570
            Value rem = (*_capacity)[e] - (*_flow)[e];
571 571
            if (!_tolerance.positive(rem)) continue;
572 572
            Node v = _graph.target(e);
573 573
            if ((*_level)[v] < level) {
574 574
              if (!_level->active(v) && v != _target) {
575 575
                _level->activate(v);
576 576
              }
577 577
              if (!_tolerance.less(rem, excess)) {
578 578
                _flow->set(e, (*_flow)[e] + excess);
579 579
                _excess->set(v, (*_excess)[v] + excess);
580 580
                excess = 0;
581 581
                goto no_more_push_1;
582 582
              } else {
583 583
                excess -= rem;
584 584
                _excess->set(v, (*_excess)[v] + rem);
585 585
                _flow->set(e, (*_capacity)[e]);
586 586
              }
587 587
            } else if (new_level > (*_level)[v]) {
588 588
              new_level = (*_level)[v];
589 589
            }
590 590
          }
591 591

	
592 592
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
593 593
            Value rem = (*_flow)[e];
594 594
            if (!_tolerance.positive(rem)) continue;
595 595
            Node v = _graph.source(e);
596 596
            if ((*_level)[v] < level) {
597 597
              if (!_level->active(v) && v != _target) {
598 598
                _level->activate(v);
599 599
              }
600 600
              if (!_tolerance.less(rem, excess)) {
601 601
                _flow->set(e, (*_flow)[e] - excess);
602 602
                _excess->set(v, (*_excess)[v] + excess);
603 603
                excess = 0;
604 604
                goto no_more_push_1;
605 605
              } else {
606 606
                excess -= rem;
607 607
                _excess->set(v, (*_excess)[v] + rem);
608 608
                _flow->set(e, 0);
609 609
              }
610 610
            } else if (new_level > (*_level)[v]) {
611 611
              new_level = (*_level)[v];
612 612
            }
613 613
          }
614 614

	
615 615
        no_more_push_1:
616 616

	
617 617
          _excess->set(n, excess);
618 618

	
619 619
          if (excess != 0) {
620 620
            if (new_level + 1 < _level->maxLevel()) {
621 621
              _level->liftHighestActive(new_level + 1);
622 622
            } else {
623 623
              _level->liftHighestActiveToTop();
624 624
            }
625 625
            if (_level->emptyLevel(level)) {
626 626
              _level->liftToTop(level);
627 627
            }
628 628
          } else {
629 629
            _level->deactivate(n);
630 630
          }
631 631

	
632 632
          n = _level->highestActive();
633 633
          level = _level->highestActiveLevel();
634 634
          --num;
635 635
        }
636 636

	
637 637
        num = _node_num * 20;
638 638
        while (num > 0 && n != INVALID) {
639 639
          Value excess = (*_excess)[n];
640 640
          int new_level = _level->maxLevel();
641 641

	
642 642
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
643 643
            Value rem = (*_capacity)[e] - (*_flow)[e];
644 644
            if (!_tolerance.positive(rem)) continue;
645 645
            Node v = _graph.target(e);
646 646
            if ((*_level)[v] < level) {
647 647
              if (!_level->active(v) && v != _target) {
648 648
                _level->activate(v);
649 649
              }
650 650
              if (!_tolerance.less(rem, excess)) {
651 651
                _flow->set(e, (*_flow)[e] + excess);
652 652
                _excess->set(v, (*_excess)[v] + excess);
653 653
                excess = 0;
654 654
                goto no_more_push_2;
655 655
              } else {
656 656
                excess -= rem;
657 657
                _excess->set(v, (*_excess)[v] + rem);
658 658
                _flow->set(e, (*_capacity)[e]);
659 659
              }
660 660
            } else if (new_level > (*_level)[v]) {
661 661
              new_level = (*_level)[v];
662 662
            }
663 663
          }
664 664

	
665 665
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
666 666
            Value rem = (*_flow)[e];
667 667
            if (!_tolerance.positive(rem)) continue;
668 668
            Node v = _graph.source(e);
669 669
            if ((*_level)[v] < level) {
670 670
              if (!_level->active(v) && v != _target) {
671 671
                _level->activate(v);
672 672
              }
673 673
              if (!_tolerance.less(rem, excess)) {
674 674
                _flow->set(e, (*_flow)[e] - excess);
675 675
                _excess->set(v, (*_excess)[v] + excess);
676 676
                excess = 0;
677 677
                goto no_more_push_2;
678 678
              } else {
679 679
                excess -= rem;
680 680
                _excess->set(v, (*_excess)[v] + rem);
681 681
                _flow->set(e, 0);
682 682
              }
683 683
            } else if (new_level > (*_level)[v]) {
684 684
              new_level = (*_level)[v];
685 685
            }
686 686
          }
687 687

	
688 688
        no_more_push_2:
689 689

	
690 690
          _excess->set(n, excess);
691 691

	
692 692
          if (excess != 0) {
693 693
            if (new_level + 1 < _level->maxLevel()) {
694 694
              _level->liftActiveOn(level, new_level + 1);
695 695
            } else {
696 696
              _level->liftActiveToTop(level);
697 697
            }
698 698
            if (_level->emptyLevel(level)) {
699 699
              _level->liftToTop(level);
700 700
            }
701 701
          } else {
702 702
            _level->deactivate(n);
703 703
          }
704 704

	
705 705
          while (level >= 0 && _level->activeFree(level)) {
706 706
            --level;
707 707
          }
708 708
          if (level == -1) {
709 709
            n = _level->highestActive();
710 710
            level = _level->highestActiveLevel();
711 711
          } else {
712 712
            n = _level->activeOn(level);
713 713
          }
714 714
          --num;
715 715
        }
716 716
      }
717 717
    }
718 718

	
719 719
    /// \brief Starts the second phase of the preflow algorithm.
720 720
    ///
721 721
    /// The preflow algorithm consists of two phases, this method runs
722 722
    /// the second phase. After calling one of the \ref init() functions
723 723
    /// and \ref startFirstPhase() and then \ref startSecondPhase(),
724 724
    /// \ref flowMap() returns a maximum flow, \ref flowValue() returns the
725 725
    /// value of a maximum flow, \ref minCut() returns a minimum cut
726 726
    /// \pre One of the \ref init() functions and \ref startFirstPhase()
727 727
    /// must be called before using this function.
728 728
    void startSecondPhase() {
729 729
      _phase = false;
730 730

	
731 731
      typename Digraph::template NodeMap<bool> reached(_graph);
732 732
      for (NodeIt n(_graph); n != INVALID; ++n) {
733 733
        reached.set(n, (*_level)[n] < _level->maxLevel());
734 734
      }
735 735

	
736 736
      _level->initStart();
737 737
      _level->initAddItem(_source);
738 738

	
739 739
      std::vector<Node> queue;
740 740
      queue.push_back(_source);
741 741
      reached.set(_source, true);
742 742

	
743 743
      while (!queue.empty()) {
744 744
        _level->initNewLevel();
745 745
        std::vector<Node> nqueue;
746 746
        for (int i = 0; i < int(queue.size()); ++i) {
747 747
          Node n = queue[i];
748 748
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
749 749
            Node v = _graph.target(e);
750 750
            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
751 751
              reached.set(v, true);
752 752
              _level->initAddItem(v);
753 753
              nqueue.push_back(v);
754 754
            }
755 755
          }
756 756
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
757 757
            Node u = _graph.source(e);
758 758
            if (!reached[u] &&
759 759
                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
760 760
              reached.set(u, true);
761 761
              _level->initAddItem(u);
762 762
              nqueue.push_back(u);
763 763
            }
764 764
          }
765 765
        }
766 766
        queue.swap(nqueue);
767 767
      }
768 768
      _level->initFinish();
769 769

	
770 770
      for (NodeIt n(_graph); n != INVALID; ++n) {
771 771
        if (!reached[n]) {
772 772
          _level->dirtyTopButOne(n);
773 773
        } else if ((*_excess)[n] > 0 && _target != n) {
774 774
          _level->activate(n);
775 775
        }
776 776
      }
777 777

	
778 778
      Node n;
779 779
      while ((n = _level->highestActive()) != INVALID) {
780 780
        Value excess = (*_excess)[n];
781 781
        int level = _level->highestActiveLevel();
782 782
        int new_level = _level->maxLevel();
783 783

	
784 784
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
785 785
          Value rem = (*_capacity)[e] - (*_flow)[e];
786 786
          if (!_tolerance.positive(rem)) continue;
787 787
          Node v = _graph.target(e);
788 788
          if ((*_level)[v] < level) {
789 789
            if (!_level->active(v) && v != _source) {
790 790
              _level->activate(v);
791 791
            }
792 792
            if (!_tolerance.less(rem, excess)) {
793 793
              _flow->set(e, (*_flow)[e] + excess);
794 794
              _excess->set(v, (*_excess)[v] + excess);
795 795
              excess = 0;
796 796
              goto no_more_push;
797 797
            } else {
798 798
              excess -= rem;
799 799
              _excess->set(v, (*_excess)[v] + rem);
800 800
              _flow->set(e, (*_capacity)[e]);
801 801
            }
802 802
          } else if (new_level > (*_level)[v]) {
803 803
            new_level = (*_level)[v];
804 804
          }
805 805
        }
806 806

	
807 807
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
808 808
          Value rem = (*_flow)[e];
809 809
          if (!_tolerance.positive(rem)) continue;
810 810
          Node v = _graph.source(e);
811 811
          if ((*_level)[v] < level) {
812 812
            if (!_level->active(v) && v != _source) {
813 813
              _level->activate(v);
814 814
            }
815 815
            if (!_tolerance.less(rem, excess)) {
816 816
              _flow->set(e, (*_flow)[e] - excess);
817 817
              _excess->set(v, (*_excess)[v] + excess);
818 818
              excess = 0;
819 819
              goto no_more_push;
820 820
            } else {
821 821
              excess -= rem;
822 822
              _excess->set(v, (*_excess)[v] + rem);
823 823
              _flow->set(e, 0);
824 824
            }
825 825
          } else if (new_level > (*_level)[v]) {
826 826
            new_level = (*_level)[v];
827 827
          }
828 828
        }
829 829

	
830 830
      no_more_push:
831 831

	
832 832
        _excess->set(n, excess);
833 833

	
834 834
        if (excess != 0) {
835 835
          if (new_level + 1 < _level->maxLevel()) {
836 836
            _level->liftHighestActive(new_level + 1);
837 837
          } else {
838 838
            // Calculation error
839 839
            _level->liftHighestActiveToTop();
840 840
          }
841 841
          if (_level->emptyLevel(level)) {
842 842
            // Calculation error
843 843
            _level->liftToTop(level);
844 844
          }
845 845
        } else {
846 846
          _level->deactivate(n);
847 847
        }
848 848

	
849 849
      }
850 850
    }
851 851

	
852 852
    /// \brief Runs the preflow algorithm.
853 853
    ///
854 854
    /// Runs the preflow algorithm.
855 855
    /// \note pf.run() is just a shortcut of the following code.
856 856
    /// \code
857 857
    ///   pf.init();
858 858
    ///   pf.startFirstPhase();
859 859
    ///   pf.startSecondPhase();
860 860
    /// \endcode
861 861
    void run() {
862 862
      init();
863 863
      startFirstPhase();
864 864
      startSecondPhase();
865 865
    }
866 866

	
867 867
    /// \brief Runs the preflow algorithm to compute the minimum cut.
868 868
    ///
869 869
    /// Runs the preflow algorithm to compute the minimum cut.
870 870
    /// \note pf.runMinCut() is just a shortcut of the following code.
871 871
    /// \code
872 872
    ///   pf.init();
873 873
    ///   pf.startFirstPhase();
874 874
    /// \endcode
875 875
    void runMinCut() {
876 876
      init();
877 877
      startFirstPhase();
878 878
    }
879 879

	
880 880
    /// @}
881 881

	
882 882
    /// \name Query Functions
883 883
    /// The results of the preflow algorithm can be obtained using these
884 884
    /// functions.\n
885 885
    /// Either one of the \ref run() "run*()" functions or one of the
886 886
    /// \ref startFirstPhase() "start*()" functions should be called
887 887
    /// before using them.
888 888

	
889 889
    ///@{
890 890

	
891 891
    /// \brief Returns the value of the maximum flow.
892 892
    ///
893 893
    /// Returns the value of the maximum flow by returning the excess
894 894
    /// of the target node. This value equals to the value of
895 895
    /// the maximum flow already after the first phase of the algorithm.
896 896
    ///
897 897
    /// \pre Either \ref run() or \ref init() must be called before
898 898
    /// using this function.
899 899
    Value flowValue() const {
900 900
      return (*_excess)[_target];
901 901
    }
902 902

	
903 903
    /// \brief Returns the flow on the given arc.
904 904
    ///
905 905
    /// Returns the flow on the given arc. This method can
906 906
    /// be called after the second phase of the algorithm.
907 907
    ///
908 908
    /// \pre Either \ref run() or \ref init() must be called before
909 909
    /// using this function.
910 910
    Value flow(const Arc& arc) const {
911 911
      return (*_flow)[arc];
912 912
    }
913 913

	
914 914
    /// \brief Returns a const reference to the flow map.
915 915
    ///
916 916
    /// Returns a const reference to the arc map storing the found flow.
917 917
    /// This method can be called after the second phase of the algorithm.
918 918
    ///
919 919
    /// \pre Either \ref run() or \ref init() must be called before
920 920
    /// using this function.
921
    const FlowMap& flowMap() {
921
    const FlowMap& flowMap() const {
922 922
      return *_flow;
923 923
    }
924 924

	
925 925
    /// \brief Returns \c true when the node is on the source side of the
926 926
    /// minimum cut.
927 927
    ///
928 928
    /// Returns true when the node is on the source side of the found
929 929
    /// minimum cut. This method can be called both after running \ref
930 930
    /// startFirstPhase() and \ref startSecondPhase().
931 931
    ///
932 932
    /// \pre Either \ref run() or \ref init() must be called before
933 933
    /// using this function.
934 934
    bool minCut(const Node& node) const {
935 935
      return ((*_level)[node] == _level->maxLevel()) == _phase;
936 936
    }
937 937

	
938 938
    /// \brief Gives back a minimum value cut.
939 939
    ///
940 940
    /// Sets \c cutMap to the characteristic vector of a minimum value
941 941
    /// cut. \c cutMap should be a \ref concepts::WriteMap "writable"
942 942
    /// node map with \c bool (or convertible) value type.
943 943
    ///
944 944
    /// This method can be called both after running \ref startFirstPhase()
945 945
    /// and \ref startSecondPhase(). The result after the second phase
946 946
    /// could be slightly different if inexact computation is used.
947 947
    ///
948 948
    /// \note This function calls \ref minCut() for each node, so it runs in
949 949
    /// \f$O(n)\f$ time.
950 950
    ///
951 951
    /// \pre Either \ref run() or \ref init() must be called before
952 952
    /// using this function.
953 953
    template <typename CutMap>
954 954
    void minCutMap(CutMap& cutMap) const {
955 955
      for (NodeIt n(_graph); n != INVALID; ++n) {
956 956
        cutMap.set(n, minCut(n));
957 957
      }
958 958
    }
959 959

	
960 960
    /// @}
961 961
  };
962 962
}
963 963

	
964 964
#endif
0 comments (0 inline)