1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/dinitz_sleator_tarjan.h Sat Nov 17 20:58:11 2007 +0000
1.3 @@ -0,0 +1,762 @@
1.4 +/* -*- C++ -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library
1.7 + *
1.8 + * Copyright (C) 2003-2007
1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 + *
1.12 + * Permission to use, modify and distribute this software is granted
1.13 + * provided that this copyright notice appears in all copies. For
1.14 + * precise terms see the accompanying LICENSE file.
1.15 + *
1.16 + * This software is provided "AS IS" with no warranty of any kind,
1.17 + * express or implied, and with no claim as to its suitability for any
1.18 + * purpose.
1.19 + *
1.20 + */
1.21 +
1.22 +#ifndef LEMON_DINITZ_SLEATOR_TARJAN_H
1.23 +#define LEMON_DINITZ_SLEATOR_TARJAN_H
1.24 +
1.25 +/// \file
1.26 +/// \ingroup max_flow
1.27 +/// \brief Implementation the dynamic tree data structure of Sleator
1.28 +/// and Tarjan.
1.29 +
1.30 +#include <lemon/time_measure.h>
1.31 +#include <queue>
1.32 +#include <lemon/graph_utils.h>
1.33 +#include <lemon/tolerance.h>
1.34 +#include <lemon/graph_adaptor.h>
1.35 +#include <lemon/bfs.h>
1.36 +#include <lemon/edge_set.h>
1.37 +#include <lemon/dynamic_tree.h>
1.38 +
1.39 +#include <vector>
1.40 +#include <limits>
1.41 +#include <fstream>
1.42 +
1.43 +
1.44 +namespace lemon {
1.45 +
1.46 + /// \brief Default traits class of DinitzSleatorTarjan class.
1.47 + ///
1.48 + /// Default traits class of DinitzSleatorTarjan class.
1.49 + /// \param _Graph Graph type.
1.50 + /// \param _CapacityMap Type of capacity map.
1.51 + template <typename _Graph, typename _CapacityMap>
1.52 + struct DinitzSleatorTarjanDefaultTraits {
1.53 +
1.54 + /// \brief The graph type the algorithm runs on.
1.55 + typedef _Graph Graph;
1.56 +
1.57 + /// \brief The type of the map that stores the edge capacities.
1.58 + ///
1.59 + /// The type of the map that stores the edge capacities.
1.60 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
1.61 + typedef _CapacityMap CapacityMap;
1.62 +
1.63 + /// \brief The type of the length of the edges.
1.64 + typedef typename CapacityMap::Value Value;
1.65 +
1.66 + /// \brief The map type that stores the flow values.
1.67 + ///
1.68 + /// The map type that stores the flow values.
1.69 + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
1.70 + typedef typename Graph::template EdgeMap<Value> FlowMap;
1.71 +
1.72 + /// \brief Instantiates a FlowMap.
1.73 + ///
1.74 + /// This function instantiates a \ref FlowMap.
1.75 + /// \param graph The graph, to which we would like to define the flow map.
1.76 + static FlowMap* createFlowMap(const Graph& graph) {
1.77 + return new FlowMap(graph);
1.78 + }
1.79 +
1.80 + /// \brief The tolerance used by the algorithm
1.81 + ///
1.82 + /// The tolerance used by the algorithm to handle inexact computation.
1.83 + typedef Tolerance<Value> Tolerance;
1.84 +
1.85 + };
1.86 +
1.87 + /// \ingroup max_flow
1.88 + ///
1.89 + /// \brief Dinitz-Sleator-Tarjan algorithms class.
1.90 + ///
1.91 + /// This class provides an implementation of the \e
1.92 + /// Dinitz-Sleator-Tarjan \e algorithm producing a flow of maximum
1.93 + /// value in a directed graph. The DinitzSleatorTarjan algorithm is
1.94 + /// the fastest known max flow algorithms wich using blocking
1.95 + /// flow. It is an improvement of the Dinitz algorithm by using the
1.96 + /// \ref DynamicTree "dynamic tree" data structure of Sleator and
1.97 + /// Tarjan.
1.98 + ///
1.99 + /// This blocking flow algorithms builds a layered graph according
1.100 + /// to \ref Bfs "breadth-first search" distance from the target node
1.101 + /// in the reversed residual graph. The layered graph contains each
1.102 + /// residual edge which steps one level down. After that the
1.103 + /// algorithm constructs a blocking flow on the layered graph and
1.104 + /// augments the overall flow with it. The number of the levels in
1.105 + /// the layered graph is strictly increasing in each augmenting
1.106 + /// phase therefore the number of the augmentings is at most
1.107 + /// \f$n-1\f$. The length of each phase is at most
1.108 + /// \f$O(m\log(n))\f$, that the overall time complexity is
1.109 + /// \f$O(nm\log(n))\f$.
1.110 + ///
1.111 + /// \param _Graph The directed graph type the algorithm runs on.
1.112 + /// \param _CapacityMap The capacity map type.
1.113 + /// \param _Traits Traits class to set various data types used by
1.114 + /// the algorithm. The default traits class is \ref
1.115 + /// DinitzSleatorTarjanDefaultTraits. See \ref
1.116 + /// DinitzSleatorTarjanDefaultTraits for the documentation of a
1.117 + /// Dinitz-Sleator-Tarjan traits class.
1.118 + ///
1.119 + /// \author Tamas Hamori and Balazs Dezso
1.120 +#ifdef DOXYGEN
1.121 + template <typename _Graph, typename _CapacityMap, typename _Traits>
1.122 +#else
1.123 + template <typename _Graph,
1.124 + typename _CapacityMap = typename _Graph::template EdgeMap<int>,
1.125 + typename _Traits =
1.126 + DinitzSleatorTarjanDefaultTraits<_Graph, _CapacityMap> >
1.127 +#endif
1.128 + class DinitzSleatorTarjan {
1.129 + public:
1.130 +
1.131 + typedef _Traits Traits;
1.132 + typedef typename Traits::Graph Graph;
1.133 + typedef typename Traits::CapacityMap CapacityMap;
1.134 + typedef typename Traits::Value Value;
1.135 +
1.136 + typedef typename Traits::FlowMap FlowMap;
1.137 + typedef typename Traits::Tolerance Tolerance;
1.138 +
1.139 +
1.140 + private:
1.141 +
1.142 + GRAPH_TYPEDEFS(typename Graph);
1.143 +
1.144 +
1.145 + typedef typename Graph::template NodeMap<int> LevelMap;
1.146 + typedef typename Graph::template NodeMap<int> IntNodeMap;
1.147 + typedef typename Graph::template NodeMap<Edge> EdgeNodeMap;
1.148 + typedef DynamicTree<Value, IntNodeMap, Tolerance, false> DynTree;
1.149 +
1.150 + private:
1.151 +
1.152 + const Graph& _graph;
1.153 + const CapacityMap* _capacity;
1.154 +
1.155 + Node _source, _target;
1.156 +
1.157 + FlowMap* _flow;
1.158 + bool _local_flow;
1.159 +
1.160 + IntNodeMap* _level;
1.161 + EdgeNodeMap* _dt_edges;
1.162 +
1.163 + IntNodeMap* _dt_index;
1.164 + DynTree* _dt;
1.165 +
1.166 + Tolerance _tolerance;
1.167 +
1.168 + Value _flow_value;
1.169 + Value _max_value;
1.170 +
1.171 +
1.172 + public:
1.173 +
1.174 + typedef DinitzSleatorTarjan Create;
1.175 +
1.176 + ///\name Named template parameters
1.177 +
1.178 + ///@{
1.179 +
1.180 + template <typename _FlowMap>
1.181 + struct DefFlowMapTraits : public Traits {
1.182 + typedef _FlowMap FlowMap;
1.183 + static FlowMap *createFlowMap(const Graph&) {
1.184 + throw UninitializedParameter();
1.185 + }
1.186 + };
1.187 +
1.188 + /// \brief \ref named-templ-param "Named parameter" for setting
1.189 + /// FlowMap type
1.190 + ///
1.191 + /// \ref named-templ-param "Named parameter" for setting FlowMap
1.192 + /// type
1.193 + template <typename _FlowMap>
1.194 + struct DefFlowMap
1.195 + : public DinitzSleatorTarjan<Graph, CapacityMap,
1.196 + DefFlowMapTraits<_FlowMap> > {
1.197 + typedef DinitzSleatorTarjan<Graph, CapacityMap,
1.198 + DefFlowMapTraits<_FlowMap> > Create;
1.199 + };
1.200 +
1.201 + template <typename _Elevator>
1.202 + struct DefElevatorTraits : public Traits {
1.203 + typedef _Elevator Elevator;
1.204 + static Elevator *createElevator(const Graph&, int) {
1.205 + throw UninitializedParameter();
1.206 + }
1.207 + };
1.208 +
1.209 + /// @}
1.210 +
1.211 + /// \brief \ref Exception for the case when the source equals the target.
1.212 + ///
1.213 + /// \ref Exception for the case when the source equals the target.
1.214 + ///
1.215 + class InvalidArgument : public lemon::LogicError {
1.216 + public:
1.217 + virtual const char* what() const throw() {
1.218 + return "lemon::DinitzSleatorTarjan::InvalidArgument";
1.219 + }
1.220 + };
1.221 +
1.222 + /// \brief The constructor of the class.
1.223 + ///
1.224 + /// The constructor of the class.
1.225 + /// \param graph The directed graph the algorithm runs on.
1.226 + /// \param capacity The capacity of the edges.
1.227 + /// \param source The source node.
1.228 + /// \param target The target node.
1.229 + DinitzSleatorTarjan(const Graph& graph, const CapacityMap& capacity,
1.230 + Node source, Node target)
1.231 + : _graph(graph), _capacity(&capacity),
1.232 + _source(source), _target(target),
1.233 + _flow(0), _local_flow(false),
1.234 + _level(0), _dt_edges(0),
1.235 + _dt_index(0), _dt(0),
1.236 + _tolerance(), _flow_value(), _max_value()
1.237 + {
1.238 + if (_source == _target) {
1.239 + throw InvalidArgument();
1.240 + }
1.241 + }
1.242 +
1.243 + /// \brief Destrcutor.
1.244 + ///
1.245 + /// Destructor.
1.246 + ~DinitzSleatorTarjan() {
1.247 + destroyStructures();
1.248 + }
1.249 +
1.250 + /// \brief Sets the capacity map.
1.251 + ///
1.252 + /// Sets the capacity map.
1.253 + /// \return \c (*this)
1.254 + DinitzSleatorTarjan& capacityMap(const CapacityMap& map) {
1.255 + _capacity = ↦
1.256 + return *this;
1.257 + }
1.258 +
1.259 + /// \brief Sets the flow map.
1.260 + ///
1.261 + /// Sets the flow map.
1.262 + /// \return \c (*this)
1.263 + DinitzSleatorTarjan& flowMap(FlowMap& map) {
1.264 + if (_local_flow) {
1.265 + delete _flow;
1.266 + _local_flow = false;
1.267 + }
1.268 + _flow = ↦
1.269 + return *this;
1.270 + }
1.271 +
1.272 + /// \brief Returns the flow map.
1.273 + ///
1.274 + /// \return The flow map.
1.275 + const FlowMap& flowMap() {
1.276 + return *_flow;
1.277 + }
1.278 +
1.279 + /// \brief Sets the source node.
1.280 + ///
1.281 + /// Sets the source node.
1.282 + /// \return \c (*this)
1.283 + DinitzSleatorTarjan& source(const Node& node) {
1.284 + _source = node;
1.285 + return *this;
1.286 + }
1.287 +
1.288 + /// \brief Sets the target node.
1.289 + ///
1.290 + /// Sets the target node.
1.291 + /// \return \c (*this)
1.292 + DinitzSleatorTarjan& target(const Node& node) {
1.293 + _target = node;
1.294 + return *this;
1.295 + }
1.296 +
1.297 + /// \brief Sets the tolerance used by algorithm.
1.298 + ///
1.299 + /// Sets the tolerance used by algorithm.
1.300 + DinitzSleatorTarjan& tolerance(const Tolerance& tolerance) const {
1.301 + _tolerance = tolerance;
1.302 + if (_dt) {
1.303 + _dt.tolerance(_tolerance);
1.304 + }
1.305 + return *this;
1.306 + }
1.307 +
1.308 + /// \brief Returns the tolerance used by algorithm.
1.309 + ///
1.310 + /// Returns the tolerance used by algorithm.
1.311 + const Tolerance& tolerance() const {
1.312 + return tolerance;
1.313 + }
1.314 +
1.315 + private:
1.316 +
1.317 + void createStructures() {
1.318 + if (!_flow) {
1.319 + _flow = Traits::createFlowMap(_graph);
1.320 + _local_flow = true;
1.321 + }
1.322 + if (!_level) {
1.323 + _level = new LevelMap(_graph);
1.324 + }
1.325 + if (!_dt_index && !_dt) {
1.326 + _dt_index = new IntNodeMap(_graph);
1.327 + _dt = new DynTree(*_dt_index, _tolerance);
1.328 + }
1.329 + if (!_dt_edges) {
1.330 + _dt_edges = new EdgeNodeMap(_graph);
1.331 + }
1.332 + _max_value = _dt->maxValue();
1.333 + }
1.334 +
1.335 + void destroyStructures() {
1.336 + if (_local_flow) {
1.337 + delete _flow;
1.338 + }
1.339 + if (_level) {
1.340 + delete _level;
1.341 + }
1.342 + if (_dt) {
1.343 + delete _dt;
1.344 + }
1.345 + if (_dt_index) {
1.346 + delete _dt_index;
1.347 + }
1.348 + if (_dt_edges) {
1.349 + delete _dt_edges;
1.350 + }
1.351 + }
1.352 +
1.353 + bool createLayeredGraph() {
1.354 +
1.355 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.356 + _level->set(n, -2);
1.357 + }
1.358 +
1.359 + int level = 0;
1.360 +
1.361 + std::vector<Node> queue;
1.362 + queue.push_back(_target);
1.363 + _level->set(_target, level);
1.364 +
1.365 + while ((*_level)[_source] == -2 && !queue.empty()) {
1.366 + std::vector<Node> nqueue;
1.367 + ++level;
1.368 +
1.369 + for (int i = 0; i < int(queue.size()); ++i) {
1.370 + Node n = queue[i];
1.371 +
1.372 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
1.373 + Node v = _graph.target(e);
1.374 + if ((*_level)[v] != -2) continue;
1.375 + Value rem = (*_flow)[e];
1.376 + if (!_tolerance.positive(rem)) continue;
1.377 + _level->set(v, level);
1.378 + nqueue.push_back(v);
1.379 + }
1.380 +
1.381 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
1.382 + Node v = _graph.source(e);
1.383 + if ((*_level)[v] != -2) continue;
1.384 + Value rem = (*_capacity)[e] - (*_flow)[e];
1.385 + if (!_tolerance.positive(rem)) continue;
1.386 + _level->set(v, level);
1.387 + nqueue.push_back(v);
1.388 + }
1.389 +
1.390 + }
1.391 + queue.swap(nqueue);
1.392 + }
1.393 +
1.394 + return (*_level)[_source] != -2;
1.395 + }
1.396 +
1.397 + void initEdges() {
1.398 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.399 + _graph.firstOut((*_dt_edges)[n], n);
1.400 + }
1.401 + }
1.402 +
1.403 +
1.404 + void augmentPath() {
1.405 + Value rem;
1.406 + Node n = _dt->findCost(_source, rem);
1.407 + _flow_value += rem;
1.408 + _dt->addCost(_source, - rem);
1.409 +
1.410 + _dt->cut(n);
1.411 + _dt->addCost(n, _max_value);
1.412 +
1.413 + Edge e = (*_dt_edges)[n];
1.414 + if (_graph.source(e) == n) {
1.415 + _flow->set(e, (*_capacity)[e]);
1.416 +
1.417 + _graph.nextOut(e);
1.418 + if (e == INVALID) {
1.419 + _graph.firstIn(e, n);
1.420 + }
1.421 + } else {
1.422 + _flow->set(e, 0);
1.423 + _graph.nextIn(e);
1.424 + }
1.425 + _dt_edges->set(n, e);
1.426 +
1.427 + }
1.428 +
1.429 + bool advance(Node n) {
1.430 + Edge e = (*_dt_edges)[n];
1.431 + if (e == INVALID) return false;
1.432 +
1.433 + Node u;
1.434 + Value rem;
1.435 + if (_graph.source(e) == n) {
1.436 + u = _graph.target(e);
1.437 + while ((*_level)[n] != (*_level)[u] + 1 ||
1.438 + !_tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
1.439 + _graph.nextOut(e);
1.440 + if (e == INVALID) break;
1.441 + u = _graph.target(e);
1.442 + }
1.443 + if (e != INVALID) {
1.444 + rem = (*_capacity)[e] - (*_flow)[e];
1.445 + } else {
1.446 + _graph.firstIn(e, n);
1.447 + if (e == INVALID) {
1.448 + _dt_edges->set(n, INVALID);
1.449 + return false;
1.450 + }
1.451 + u = _graph.source(e);
1.452 + while ((*_level)[n] != (*_level)[u] + 1 ||
1.453 + !_tolerance.positive((*_flow)[e])) {
1.454 + _graph.nextIn(e);
1.455 + if (e == INVALID) {
1.456 + _dt_edges->set(n, INVALID);
1.457 + return false;
1.458 + }
1.459 + u = _graph.source(e);
1.460 + }
1.461 + rem = (*_flow)[e];
1.462 + }
1.463 + } else {
1.464 + u = _graph.source(e);
1.465 + while ((*_level)[n] != (*_level)[u] + 1 ||
1.466 + !_tolerance.positive((*_flow)[e])) {
1.467 + _graph.nextIn(e);
1.468 + if (e == INVALID) {
1.469 + _dt_edges->set(n, INVALID);
1.470 + return false;
1.471 + }
1.472 + u = _graph.source(e);
1.473 + }
1.474 + rem = (*_flow)[e];
1.475 + }
1.476 +
1.477 + _dt->addCost(n, - std::numeric_limits<Value>::max());
1.478 + _dt->addCost(n, rem);
1.479 + _dt->link(n, u);
1.480 + _dt_edges->set(n, e);
1.481 + return true;
1.482 + }
1.483 +
1.484 + void retreat(Node n) {
1.485 + _level->set(n, -1);
1.486 +
1.487 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
1.488 + Node u = _graph.target(e);
1.489 + if ((*_dt_edges)[u] == e && _dt->findRoot(u) == n) {
1.490 + Value rem;
1.491 + _dt->findCost(u, rem);
1.492 + _flow->set(e, rem);
1.493 + _dt->cut(u);
1.494 + _dt->addCost(u, - rem);
1.495 + _dt->addCost(u, _max_value);
1.496 + }
1.497 + }
1.498 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
1.499 + Node u = _graph.source(e);
1.500 + if ((*_dt_edges)[u] == e && _dt->findRoot(u) == n) {
1.501 + Value rem;
1.502 + _dt->findCost(u, rem);
1.503 + _flow->set(e, (*_capacity)[e] - rem);
1.504 + _dt->cut(u);
1.505 + _dt->addCost(u, - rem);
1.506 + _dt->addCost(u, _max_value);
1.507 + }
1.508 + }
1.509 + }
1.510 +
1.511 + void extractTrees() {
1.512 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.513 +
1.514 + Node w = _dt->findRoot(n);
1.515 +
1.516 + while (w != n) {
1.517 +
1.518 + Value rem;
1.519 + Node u = _dt->findCost(n, rem);
1.520 +
1.521 + _dt->cut(u);
1.522 + _dt->addCost(u, - rem);
1.523 + _dt->addCost(u, _max_value);
1.524 +
1.525 + Edge e = (*_dt_edges)[u];
1.526 + _dt_edges->set(u, INVALID);
1.527 +
1.528 + if (u == _graph.source(e)) {
1.529 + _flow->set(e, (*_capacity)[e] - rem);
1.530 + } else {
1.531 + _flow->set(e, rem);
1.532 + }
1.533 +
1.534 + w = _dt->findRoot(n);
1.535 + }
1.536 + }
1.537 + }
1.538 +
1.539 +
1.540 + public:
1.541 +
1.542 + /// \name Execution control The simplest way to execute the
1.543 + /// algorithm is to use the \c run() member functions.
1.544 + /// \n
1.545 + /// If you need more control on initial solution or
1.546 + /// execution then you have to call one \ref init() function and then
1.547 + /// the start() or multiple times the \c augment() member function.
1.548 +
1.549 + ///@{
1.550 +
1.551 + /// \brief Initializes the algorithm
1.552 + ///
1.553 + /// It sets the flow to empty flow.
1.554 + void init() {
1.555 + createStructures();
1.556 +
1.557 + _dt->clear();
1.558 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.559 + _dt->makeTree(n);
1.560 + _dt->addCost(n, _max_value);
1.561 + }
1.562 +
1.563 + for (EdgeIt it(_graph); it != INVALID; ++it) {
1.564 + _flow->set(it, 0);
1.565 + }
1.566 + _flow_value = 0;
1.567 + }
1.568 +
1.569 + /// \brief Initializes the algorithm
1.570 + ///
1.571 + /// Initializes the flow to the \c flowMap. The \c flowMap should
1.572 + /// contain a feasible flow, ie. in each node excluding the source
1.573 + /// and the target the incoming flow should be equal to the
1.574 + /// outgoing flow.
1.575 + template <typename FlowMap>
1.576 + void flowInit(const FlowMap& flowMap) {
1.577 + createStructures();
1.578 +
1.579 + _dt->clear();
1.580 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.581 + _dt->makeTree(n);
1.582 + _dt->addCost(n, _max_value);
1.583 + }
1.584 +
1.585 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.586 + _flow->set(e, flowMap[e]);
1.587 + }
1.588 + _flow_value = 0;
1.589 + for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
1.590 + _flow_value += (*_flow)[jt];
1.591 + }
1.592 + for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
1.593 + _flow_value -= (*_flow)[jt];
1.594 + }
1.595 + }
1.596 +
1.597 + /// \brief Initializes the algorithm
1.598 + ///
1.599 + /// Initializes the flow to the \c flowMap. The \c flowMap should
1.600 + /// contain a feasible flow, ie. in each node excluding the source
1.601 + /// and the target the incoming flow should be equal to the
1.602 + /// outgoing flow.
1.603 + /// \return %False when the given flowMap does not contain
1.604 + /// feasible flow.
1.605 + template <typename FlowMap>
1.606 + bool checkedFlowInit(const FlowMap& flowMap) {
1.607 + createStructures();
1.608 +
1.609 + _dt->clear();
1.610 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.611 + _dt->makeTree(n);
1.612 + _dt->addCost(n, _max_value);
1.613 + }
1.614 +
1.615 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.616 + _flow->set(e, flowMap[e]);
1.617 + }
1.618 + for (NodeIt it(_graph); it != INVALID; ++it) {
1.619 + if (it == _source || it == _target) continue;
1.620 + Value outFlow = 0;
1.621 + for (OutEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
1.622 + outFlow += (*_flow)[jt];
1.623 + }
1.624 + Value inFlow = 0;
1.625 + for (InEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
1.626 + inFlow += (*_flow)[jt];
1.627 + }
1.628 + if (_tolerance.different(outFlow, inFlow)) {
1.629 + return false;
1.630 + }
1.631 + }
1.632 + for (EdgeIt it(_graph); it != INVALID; ++it) {
1.633 + if (_tolerance.less((*_flow)[it], 0)) return false;
1.634 + if (_tolerance.less((*_capacity)[it], (*_flow)[it])) return false;
1.635 + }
1.636 + _flow_value = 0;
1.637 + for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
1.638 + _flow_value += (*_flow)[jt];
1.639 + }
1.640 + for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
1.641 + _flow_value -= (*_flow)[jt];
1.642 + }
1.643 + return true;
1.644 + }
1.645 +
1.646 + /// \brief Executes the algorithm
1.647 + ///
1.648 + /// It runs augmenting phases by adding blocking flow until the
1.649 + /// optimal solution is reached.
1.650 + void start() {
1.651 + while (augment());
1.652 + }
1.653 +
1.654 + /// \brief Augments the flow with a blocking flow on a layered
1.655 + /// graph.
1.656 + ///
1.657 + /// This function builds a layered graph and then find a blocking
1.658 + /// flow on this graph. The number of the levels in the layered
1.659 + /// graph is strictly increasing in each augmenting phase
1.660 + /// therefore the number of the augmentings is at most \f$ n-1
1.661 + /// \f$. The length of each phase is at most \f$ O(m \log(n))
1.662 + /// \f$, that the overall time complexity is \f$ O(nm \log(n)) \f$.
1.663 + /// \return %False when there is not residual path between the
1.664 + /// source and the target so the current flow is a feasible and
1.665 + /// optimal solution.
1.666 + bool augment() {
1.667 + Node n;
1.668 +
1.669 + if (createLayeredGraph()) {
1.670 +
1.671 + Timer bf_timer;
1.672 + initEdges();
1.673 +
1.674 + n = _dt->findRoot(_source);
1.675 + while (true) {
1.676 + Edge e;
1.677 + if (n == _target) {
1.678 + augmentPath();
1.679 + } else if (!advance(n)) {
1.680 + if (n != _source) {
1.681 + retreat(n);
1.682 + } else {
1.683 + break;
1.684 + }
1.685 + }
1.686 + n = _dt->findRoot(_source);
1.687 + }
1.688 + extractTrees();
1.689 +
1.690 + return true;
1.691 + } else {
1.692 + return false;
1.693 + }
1.694 + }
1.695 +
1.696 + /// \brief runs the algorithm.
1.697 + ///
1.698 + /// It is just a shorthand for:
1.699 + ///
1.700 + ///\code
1.701 + /// ek.init();
1.702 + /// ek.start();
1.703 + ///\endcode
1.704 + void run() {
1.705 + init();
1.706 + start();
1.707 + }
1.708 +
1.709 + /// @}
1.710 +
1.711 + /// \name Query Functions
1.712 + /// The result of the %Dijkstra algorithm can be obtained using these
1.713 + /// functions.\n
1.714 + /// Before the use of these functions,
1.715 + /// either run() or start() must be called.
1.716 +
1.717 + ///@{
1.718 +
1.719 + /// \brief Returns the value of the maximum flow.
1.720 + ///
1.721 + /// Returns the value of the maximum flow by returning the excess
1.722 + /// of the target node \c t. This value equals to the value of
1.723 + /// the maximum flow already after the first phase.
1.724 + Value flowValue() const {
1.725 + return _flow_value;
1.726 + }
1.727 +
1.728 +
1.729 + /// \brief Returns the flow on the edge.
1.730 + ///
1.731 + /// Sets the \c flowMap to the flow on the edges. This method can
1.732 + /// be called after the second phase of algorithm.
1.733 + Value flow(const Edge& edge) const {
1.734 + return (*_flow)[edge];
1.735 + }
1.736 +
1.737 + /// \brief Returns true when the node is on the source side of minimum cut.
1.738 + ///
1.739 +
1.740 + /// Returns true when the node is on the source side of minimum
1.741 + /// cut. This method can be called both after running \ref
1.742 + /// startFirstPhase() and \ref startSecondPhase().
1.743 + bool minCut(const Node& node) const {
1.744 + return (*_level)[node] == -2;
1.745 + }
1.746 +
1.747 + /// \brief Returns a minimum value cut.
1.748 + ///
1.749 + /// Sets \c cut to the characteristic vector of a minimum value cut
1.750 + /// It simply calls the minMinCut member.
1.751 + /// \retval cut Write node bool map.
1.752 + template <typename CutMap>
1.753 + void minCutMap(CutMap& cutMap) const {
1.754 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.755 + cutMap.set(n, (*_level)[n] == -2);
1.756 + }
1.757 + cutMap.set(_source, true);
1.758 + }
1.759 +
1.760 + /// @}
1.761 +
1.762 + };
1.763 +}
1.764 +
1.765 +#endif