[Lemon-commits] deba: r3380 - in lemon/trunk: benchmark demo doc doc/images lemon test
Lemon SVN
svn at lemon.cs.elte.hu
Sat Nov 17 21:58:20 CET 2007
Author: deba
Date: Sat Nov 17 21:58:11 2007
New Revision: 3380
Added:
lemon/trunk/lemon/dinitz_sleator_tarjan.h
lemon/trunk/lemon/dynamic_tree.h
lemon/trunk/lemon/goldberg_tarjan.h
Removed:
lemon/trunk/doc/images/flow.eps
lemon/trunk/doc/images/flow.png
Modified:
lemon/trunk/benchmark/hcube.cc
lemon/trunk/demo/disjoint_paths_demo.cc
lemon/trunk/demo/sub_graph_adaptor_demo.cc
lemon/trunk/doc/groups.dox
lemon/trunk/lemon/Makefile.am
lemon/trunk/lemon/edmonds_karp.h
lemon/trunk/lemon/preflow.h
lemon/trunk/test/preflow_test.cc
Log:
Redesign the maximum flow algorithms
Redesigned interface
Preflow changed to use elevator
Edmonds-Karp does not use the ResGraphAdaptor
Goldberg-Tarjan algorithm (Preflow with Dynamic Trees)
Dinitz-Sleator-Tarjan (Blocking flow with Dynamic Tree)
Modified: lemon/trunk/benchmark/hcube.cc
==============================================================================
--- lemon/trunk/benchmark/hcube.cc (original)
+++ lemon/trunk/benchmark/hcube.cc Sat Nov 17 21:58:11 2007
@@ -112,9 +112,9 @@
{
Graph::EdgeMap<int> flow(G);
- Preflow<Graph,int> MF(G,nodes[0],nodes[1<<dim-1],map,flow);
+ Preflow<Graph> MF(G,map,nodes[0],nodes[1<<dim-1]);
for(int i=0;i<10;i++)
- MF.run(MF.NO_FLOW);
+ MF.run();
}
PrintTime("PREFLOW",T);
Modified: lemon/trunk/demo/disjoint_paths_demo.cc
==============================================================================
--- lemon/trunk/demo/disjoint_paths_demo.cc (original)
+++ lemon/trunk/demo/disjoint_paths_demo.cc Sat Nov 17 21:58:11 2007
@@ -63,7 +63,7 @@
Graph::EdgeMap<int> flow(graph);
- Preflow<Graph, int, Capacity> preflow(graph, source, target, capacity, flow);
+ Preflow<Graph, Capacity> preflow(graph, capacity, source, target);
preflow.run();
@@ -72,7 +72,7 @@
graphToEps(graph, "edge_disjoint_paths.eps").
title("edge disjoint paths").scaleToA4().
copyright("(C) 2003-2007 LEMON Project").drawArrows().
- edgeColors(composeMap(functorMap(color), flow)).
+ edgeColors(composeMap(functorMap(color), preflow.flowMap())).
coords(coords).run();
@@ -87,9 +87,9 @@
SGraph::EdgeMap<int> sflow(sgraph);
- Preflow<SGraph, int, SCapacity> spreflow(sgraph, SGraph::outNode(source),
- SGraph::inNode(target),
- scapacity, sflow);
+ Preflow<SGraph, SCapacity> spreflow(sgraph, scapacity,
+ SGraph::outNode(source),
+ SGraph::inNode(target));
spreflow.run();
@@ -99,7 +99,7 @@
graphToEps(sgraph, "node_disjoint_paths.eps").
title("node disjoint paths").scaleToA4().
copyright("(C) 2003-2007 LEMON Project").drawArrows().
- edgeColors(composeMap(functorMap(color), sflow)).
+ edgeColors(composeMap(functorMap(color), spreflow.flowMap())).
coords(SGraph::combinedNodeMap(coords,
shiftMap(coords,
dim2::Point<double>(5, 0)))).
Modified: lemon/trunk/demo/sub_graph_adaptor_demo.cc
==============================================================================
--- lemon/trunk/demo/sub_graph_adaptor_demo.cc (original)
+++ lemon/trunk/demo/sub_graph_adaptor_demo.cc Sat Nov 17 21:58:11 2007
@@ -109,10 +109,8 @@
SubGW gw(g, tight_edge_filter);
ConstMap<Edge, int> const_1_map(1);
- Graph::EdgeMap<int> flow(g, 0);
// Max flow between s and t in the graph of tight edges.
- Preflow<SubGW, int, ConstMap<Edge, int>, Graph::EdgeMap<int> >
- preflow(gw, s, t, const_1_map, flow);
+ Preflow<SubGW, ConstMap<Edge, int> > preflow(gw, const_1_map, s, t);
preflow.run();
cout << "maximum number of edge-disjoint shortest paths: "
@@ -120,7 +118,7 @@
cout << "edges of the maximum number of edge-disjoint shortest s-t paths: "
<< endl;
for(EdgeIt e(g); e!=INVALID; ++e)
- if (flow[e])
+ if (preflow.flow(e) != 0)
cout << " " << g.id(e) << ", "
<< g.id(g.source(e)) << "--"
<< length[e] << "->" << g.id(g.target(e)) << endl;
Modified: lemon/trunk/doc/groups.dox
==============================================================================
--- lemon/trunk/doc/groups.dox (original)
+++ lemon/trunk/doc/groups.dox Sat Nov 17 21:58:11 2007
@@ -223,8 +223,27 @@
This group describes the algorithms for finding maximum flows and
feasible circulations.
-\image html flow.png
-\image latex flow.eps "Graph flow" width=\textwidth
+The maximum flow problem is to find a flow between a single-source and
+single-target that is maximum. Formally, there is \f$G=(V,A)\f$
+directed graph, an \f$c_a:A\rightarrow\mathbf{R}^+_0\f$ capacity
+function and given \f$s, t \in V\f$ source and target node. The
+maximum flow is the solution of the next optimization problem:
+
+\f[ 0 \le f_a \le c_a \f]
+\f[ \sum_{v\in\delta^{-}(u)}f_{vu}=\sum_{v\in\delta^{+}(u)}f_{uv} \quad u \in V \setminus \{s,t\}\f]
+\f[ \max \sum_{v\in\delta^{+}(s)}f_{uv} - \sum_{v\in\delta^{-}(s)}f_{vu}\f]
+
+The lemon contains several algorithms for solve maximum flow problems:
+- \ref lemon::EdmondsKarp "Edmonds-Karp"
+- \ref lemon::Preflow "Goldberg's Preflow algorithm"
+- \ref lemon::DinitzSleatorTarjan "Dinitz's blocking flow algorithm with dynamic tree"
+- \ref lemon::GoldbergTarjan "Preflow algorithm with dynamic trees"
+
+In most cases the \ref lemon::Preflow "preflow" algorithm provides the
+fastest method to compute the maximum flow. All impelementations
+provides functions for query the minimum cut, which is the dual linear
+programming probelm of the maximum flow.
+
*/
/**
Modified: lemon/trunk/lemon/Makefile.am
==============================================================================
--- lemon/trunk/lemon/Makefile.am (original)
+++ lemon/trunk/lemon/Makefile.am Sat Nov 17 21:58:11 2007
@@ -52,9 +52,11 @@
lemon/dag_shortest_path.h \
lemon/dfs.h \
lemon/dijkstra.h \
+ lemon/dinitz_sleator_tarjan.h \
lemon/dist_log.h \
lemon/dim2.h \
lemon/dimacs.h \
+ lemon/dynamic_tree.h \
lemon/edge_set.h \
lemon/edmonds_karp.h \
lemon/elevator.h \
@@ -71,6 +73,7 @@
lemon/graph_utils.h \
lemon/graph_writer.h \
lemon/grid_ugraph.h \
+ lemon/goldberg_tarjan.h \
lemon/hao_orlin.h \
lemon/hypercube_graph.h \
lemon/iterable_maps.h \
Added: lemon/trunk/lemon/dinitz_sleator_tarjan.h
==============================================================================
--- (empty file)
+++ lemon/trunk/lemon/dinitz_sleator_tarjan.h Sat Nov 17 21:58:11 2007
@@ -0,0 +1,762 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_DINITZ_SLEATOR_TARJAN_H
+#define LEMON_DINITZ_SLEATOR_TARJAN_H
+
+/// \file
+/// \ingroup max_flow
+/// \brief Implementation the dynamic tree data structure of Sleator
+/// and Tarjan.
+
+#include <lemon/time_measure.h>
+#include <queue>
+#include <lemon/graph_utils.h>
+#include <lemon/tolerance.h>
+#include <lemon/graph_adaptor.h>
+#include <lemon/bfs.h>
+#include <lemon/edge_set.h>
+#include <lemon/dynamic_tree.h>
+
+#include <vector>
+#include <limits>
+#include <fstream>
+
+
+namespace lemon {
+
+ /// \brief Default traits class of DinitzSleatorTarjan class.
+ ///
+ /// Default traits class of DinitzSleatorTarjan class.
+ /// \param _Graph Graph type.
+ /// \param _CapacityMap Type of capacity map.
+ template <typename _Graph, typename _CapacityMap>
+ struct DinitzSleatorTarjanDefaultTraits {
+
+ /// \brief The graph type the algorithm runs on.
+ typedef _Graph Graph;
+
+ /// \brief The type of the map that stores the edge capacities.
+ ///
+ /// The type of the map that stores the edge capacities.
+ /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
+ typedef _CapacityMap CapacityMap;
+
+ /// \brief The type of the length of the edges.
+ typedef typename CapacityMap::Value Value;
+
+ /// \brief The map type that stores the flow values.
+ ///
+ /// The map type that stores the flow values.
+ /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
+ typedef typename Graph::template EdgeMap<Value> FlowMap;
+
+ /// \brief Instantiates a FlowMap.
+ ///
+ /// This function instantiates a \ref FlowMap.
+ /// \param graph The graph, to which we would like to define the flow map.
+ static FlowMap* createFlowMap(const Graph& graph) {
+ return new FlowMap(graph);
+ }
+
+ /// \brief The tolerance used by the algorithm
+ ///
+ /// The tolerance used by the algorithm to handle inexact computation.
+ typedef Tolerance<Value> Tolerance;
+
+ };
+
+ /// \ingroup max_flow
+ ///
+ /// \brief Dinitz-Sleator-Tarjan algorithms class.
+ ///
+ /// This class provides an implementation of the \e
+ /// Dinitz-Sleator-Tarjan \e algorithm producing a flow of maximum
+ /// value in a directed graph. The DinitzSleatorTarjan algorithm is
+ /// the fastest known max flow algorithms wich using blocking
+ /// flow. It is an improvement of the Dinitz algorithm by using the
+ /// \ref DynamicTree "dynamic tree" data structure of Sleator and
+ /// Tarjan.
+ ///
+ /// This blocking flow algorithms builds a layered graph according
+ /// to \ref Bfs "breadth-first search" distance from the target node
+ /// in the reversed residual graph. The layered graph contains each
+ /// residual edge which steps one level down. After that the
+ /// algorithm constructs a blocking flow on the layered graph and
+ /// augments the overall flow with it. The number of the levels in
+ /// the layered graph is strictly increasing in each augmenting
+ /// phase therefore the number of the augmentings is at most
+ /// \f$n-1\f$. The length of each phase is at most
+ /// \f$O(m\log(n))\f$, that the overall time complexity is
+ /// \f$O(nm\log(n))\f$.
+ ///
+ /// \param _Graph The directed graph type the algorithm runs on.
+ /// \param _CapacityMap The capacity map type.
+ /// \param _Traits Traits class to set various data types used by
+ /// the algorithm. The default traits class is \ref
+ /// DinitzSleatorTarjanDefaultTraits. See \ref
+ /// DinitzSleatorTarjanDefaultTraits for the documentation of a
+ /// Dinitz-Sleator-Tarjan traits class.
+ ///
+ /// \author Tamas Hamori and Balazs Dezso
+#ifdef DOXYGEN
+ template <typename _Graph, typename _CapacityMap, typename _Traits>
+#else
+ template <typename _Graph,
+ typename _CapacityMap = typename _Graph::template EdgeMap<int>,
+ typename _Traits =
+ DinitzSleatorTarjanDefaultTraits<_Graph, _CapacityMap> >
+#endif
+ class DinitzSleatorTarjan {
+ public:
+
+ typedef _Traits Traits;
+ typedef typename Traits::Graph Graph;
+ typedef typename Traits::CapacityMap CapacityMap;
+ typedef typename Traits::Value Value;
+
+ typedef typename Traits::FlowMap FlowMap;
+ typedef typename Traits::Tolerance Tolerance;
+
+
+ private:
+
+ GRAPH_TYPEDEFS(typename Graph);
+
+
+ typedef typename Graph::template NodeMap<int> LevelMap;
+ typedef typename Graph::template NodeMap<int> IntNodeMap;
+ typedef typename Graph::template NodeMap<Edge> EdgeNodeMap;
+ typedef DynamicTree<Value, IntNodeMap, Tolerance, false> DynTree;
+
+ private:
+
+ const Graph& _graph;
+ const CapacityMap* _capacity;
+
+ Node _source, _target;
+
+ FlowMap* _flow;
+ bool _local_flow;
+
+ IntNodeMap* _level;
+ EdgeNodeMap* _dt_edges;
+
+ IntNodeMap* _dt_index;
+ DynTree* _dt;
+
+ Tolerance _tolerance;
+
+ Value _flow_value;
+ Value _max_value;
+
+
+ public:
+
+ typedef DinitzSleatorTarjan Create;
+
+ ///\name Named template parameters
+
+ ///@{
+
+ template <typename _FlowMap>
+ struct DefFlowMapTraits : public Traits {
+ typedef _FlowMap FlowMap;
+ static FlowMap *createFlowMap(const Graph&) {
+ throw UninitializedParameter();
+ }
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting
+ /// FlowMap type
+ ///
+ /// \ref named-templ-param "Named parameter" for setting FlowMap
+ /// type
+ template <typename _FlowMap>
+ struct DefFlowMap
+ : public DinitzSleatorTarjan<Graph, CapacityMap,
+ DefFlowMapTraits<_FlowMap> > {
+ typedef DinitzSleatorTarjan<Graph, CapacityMap,
+ DefFlowMapTraits<_FlowMap> > Create;
+ };
+
+ template <typename _Elevator>
+ struct DefElevatorTraits : public Traits {
+ typedef _Elevator Elevator;
+ static Elevator *createElevator(const Graph&, int) {
+ throw UninitializedParameter();
+ }
+ };
+
+ /// @}
+
+ /// \brief \ref Exception for the case when the source equals the target.
+ ///
+ /// \ref Exception for the case when the source equals the target.
+ ///
+ class InvalidArgument : public lemon::LogicError {
+ public:
+ virtual const char* what() const throw() {
+ return "lemon::DinitzSleatorTarjan::InvalidArgument";
+ }
+ };
+
+ /// \brief The constructor of the class.
+ ///
+ /// The constructor of the class.
+ /// \param graph The directed graph the algorithm runs on.
+ /// \param capacity The capacity of the edges.
+ /// \param source The source node.
+ /// \param target The target node.
+ DinitzSleatorTarjan(const Graph& graph, const CapacityMap& capacity,
+ Node source, Node target)
+ : _graph(graph), _capacity(&capacity),
+ _source(source), _target(target),
+ _flow(0), _local_flow(false),
+ _level(0), _dt_edges(0),
+ _dt_index(0), _dt(0),
+ _tolerance(), _flow_value(), _max_value()
+ {
+ if (_source == _target) {
+ throw InvalidArgument();
+ }
+ }
+
+ /// \brief Destrcutor.
+ ///
+ /// Destructor.
+ ~DinitzSleatorTarjan() {
+ destroyStructures();
+ }
+
+ /// \brief Sets the capacity map.
+ ///
+ /// Sets the capacity map.
+ /// \return \c (*this)
+ DinitzSleatorTarjan& capacityMap(const CapacityMap& map) {
+ _capacity = ↦
+ return *this;
+ }
+
+ /// \brief Sets the flow map.
+ ///
+ /// Sets the flow map.
+ /// \return \c (*this)
+ DinitzSleatorTarjan& flowMap(FlowMap& map) {
+ if (_local_flow) {
+ delete _flow;
+ _local_flow = false;
+ }
+ _flow = ↦
+ return *this;
+ }
+
+ /// \brief Returns the flow map.
+ ///
+ /// \return The flow map.
+ const FlowMap& flowMap() {
+ return *_flow;
+ }
+
+ /// \brief Sets the source node.
+ ///
+ /// Sets the source node.
+ /// \return \c (*this)
+ DinitzSleatorTarjan& source(const Node& node) {
+ _source = node;
+ return *this;
+ }
+
+ /// \brief Sets the target node.
+ ///
+ /// Sets the target node.
+ /// \return \c (*this)
+ DinitzSleatorTarjan& target(const Node& node) {
+ _target = node;
+ return *this;
+ }
+
+ /// \brief Sets the tolerance used by algorithm.
+ ///
+ /// Sets the tolerance used by algorithm.
+ DinitzSleatorTarjan& tolerance(const Tolerance& tolerance) const {
+ _tolerance = tolerance;
+ if (_dt) {
+ _dt.tolerance(_tolerance);
+ }
+ return *this;
+ }
+
+ /// \brief Returns the tolerance used by algorithm.
+ ///
+ /// Returns the tolerance used by algorithm.
+ const Tolerance& tolerance() const {
+ return tolerance;
+ }
+
+ private:
+
+ void createStructures() {
+ if (!_flow) {
+ _flow = Traits::createFlowMap(_graph);
+ _local_flow = true;
+ }
+ if (!_level) {
+ _level = new LevelMap(_graph);
+ }
+ if (!_dt_index && !_dt) {
+ _dt_index = new IntNodeMap(_graph);
+ _dt = new DynTree(*_dt_index, _tolerance);
+ }
+ if (!_dt_edges) {
+ _dt_edges = new EdgeNodeMap(_graph);
+ }
+ _max_value = _dt->maxValue();
+ }
+
+ void destroyStructures() {
+ if (_local_flow) {
+ delete _flow;
+ }
+ if (_level) {
+ delete _level;
+ }
+ if (_dt) {
+ delete _dt;
+ }
+ if (_dt_index) {
+ delete _dt_index;
+ }
+ if (_dt_edges) {
+ delete _dt_edges;
+ }
+ }
+
+ bool createLayeredGraph() {
+
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ _level->set(n, -2);
+ }
+
+ int level = 0;
+
+ std::vector<Node> queue;
+ queue.push_back(_target);
+ _level->set(_target, level);
+
+ while ((*_level)[_source] == -2 && !queue.empty()) {
+ std::vector<Node> nqueue;
+ ++level;
+
+ for (int i = 0; i < int(queue.size()); ++i) {
+ Node n = queue[i];
+
+ for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Node v = _graph.target(e);
+ if ((*_level)[v] != -2) continue;
+ Value rem = (*_flow)[e];
+ if (!_tolerance.positive(rem)) continue;
+ _level->set(v, level);
+ nqueue.push_back(v);
+ }
+
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Node v = _graph.source(e);
+ if ((*_level)[v] != -2) continue;
+ Value rem = (*_capacity)[e] - (*_flow)[e];
+ if (!_tolerance.positive(rem)) continue;
+ _level->set(v, level);
+ nqueue.push_back(v);
+ }
+
+ }
+ queue.swap(nqueue);
+ }
+
+ return (*_level)[_source] != -2;
+ }
+
+ void initEdges() {
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ _graph.firstOut((*_dt_edges)[n], n);
+ }
+ }
+
+
+ void augmentPath() {
+ Value rem;
+ Node n = _dt->findCost(_source, rem);
+ _flow_value += rem;
+ _dt->addCost(_source, - rem);
+
+ _dt->cut(n);
+ _dt->addCost(n, _max_value);
+
+ Edge e = (*_dt_edges)[n];
+ if (_graph.source(e) == n) {
+ _flow->set(e, (*_capacity)[e]);
+
+ _graph.nextOut(e);
+ if (e == INVALID) {
+ _graph.firstIn(e, n);
+ }
+ } else {
+ _flow->set(e, 0);
+ _graph.nextIn(e);
+ }
+ _dt_edges->set(n, e);
+
+ }
+
+ bool advance(Node n) {
+ Edge e = (*_dt_edges)[n];
+ if (e == INVALID) return false;
+
+ Node u;
+ Value rem;
+ if (_graph.source(e) == n) {
+ u = _graph.target(e);
+ while ((*_level)[n] != (*_level)[u] + 1 ||
+ !_tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
+ _graph.nextOut(e);
+ if (e == INVALID) break;
+ u = _graph.target(e);
+ }
+ if (e != INVALID) {
+ rem = (*_capacity)[e] - (*_flow)[e];
+ } else {
+ _graph.firstIn(e, n);
+ if (e == INVALID) {
+ _dt_edges->set(n, INVALID);
+ return false;
+ }
+ u = _graph.source(e);
+ while ((*_level)[n] != (*_level)[u] + 1 ||
+ !_tolerance.positive((*_flow)[e])) {
+ _graph.nextIn(e);
+ if (e == INVALID) {
+ _dt_edges->set(n, INVALID);
+ return false;
+ }
+ u = _graph.source(e);
+ }
+ rem = (*_flow)[e];
+ }
+ } else {
+ u = _graph.source(e);
+ while ((*_level)[n] != (*_level)[u] + 1 ||
+ !_tolerance.positive((*_flow)[e])) {
+ _graph.nextIn(e);
+ if (e == INVALID) {
+ _dt_edges->set(n, INVALID);
+ return false;
+ }
+ u = _graph.source(e);
+ }
+ rem = (*_flow)[e];
+ }
+
+ _dt->addCost(n, - std::numeric_limits<Value>::max());
+ _dt->addCost(n, rem);
+ _dt->link(n, u);
+ _dt_edges->set(n, e);
+ return true;
+ }
+
+ void retreat(Node n) {
+ _level->set(n, -1);
+
+ for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Node u = _graph.target(e);
+ if ((*_dt_edges)[u] == e && _dt->findRoot(u) == n) {
+ Value rem;
+ _dt->findCost(u, rem);
+ _flow->set(e, rem);
+ _dt->cut(u);
+ _dt->addCost(u, - rem);
+ _dt->addCost(u, _max_value);
+ }
+ }
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Node u = _graph.source(e);
+ if ((*_dt_edges)[u] == e && _dt->findRoot(u) == n) {
+ Value rem;
+ _dt->findCost(u, rem);
+ _flow->set(e, (*_capacity)[e] - rem);
+ _dt->cut(u);
+ _dt->addCost(u, - rem);
+ _dt->addCost(u, _max_value);
+ }
+ }
+ }
+
+ void extractTrees() {
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+
+ Node w = _dt->findRoot(n);
+
+ while (w != n) {
+
+ Value rem;
+ Node u = _dt->findCost(n, rem);
+
+ _dt->cut(u);
+ _dt->addCost(u, - rem);
+ _dt->addCost(u, _max_value);
+
+ Edge e = (*_dt_edges)[u];
+ _dt_edges->set(u, INVALID);
+
+ if (u == _graph.source(e)) {
+ _flow->set(e, (*_capacity)[e] - rem);
+ } else {
+ _flow->set(e, rem);
+ }
+
+ w = _dt->findRoot(n);
+ }
+ }
+ }
+
+
+ public:
+
+ /// \name Execution control The simplest way to execute the
+ /// algorithm is to use the \c run() member functions.
+ /// \n
+ /// If you need more control on initial solution or
+ /// execution then you have to call one \ref init() function and then
+ /// the start() or multiple times the \c augment() member function.
+
+ ///@{
+
+ /// \brief Initializes the algorithm
+ ///
+ /// It sets the flow to empty flow.
+ void init() {
+ createStructures();
+
+ _dt->clear();
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ _dt->makeTree(n);
+ _dt->addCost(n, _max_value);
+ }
+
+ for (EdgeIt it(_graph); it != INVALID; ++it) {
+ _flow->set(it, 0);
+ }
+ _flow_value = 0;
+ }
+
+ /// \brief Initializes the algorithm
+ ///
+ /// Initializes the flow to the \c flowMap. The \c flowMap should
+ /// contain a feasible flow, ie. in each node excluding the source
+ /// and the target the incoming flow should be equal to the
+ /// outgoing flow.
+ template <typename FlowMap>
+ void flowInit(const FlowMap& flowMap) {
+ createStructures();
+
+ _dt->clear();
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ _dt->makeTree(n);
+ _dt->addCost(n, _max_value);
+ }
+
+ for (EdgeIt e(_graph); e != INVALID; ++e) {
+ _flow->set(e, flowMap[e]);
+ }
+ _flow_value = 0;
+ for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
+ _flow_value += (*_flow)[jt];
+ }
+ for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
+ _flow_value -= (*_flow)[jt];
+ }
+ }
+
+ /// \brief Initializes the algorithm
+ ///
+ /// Initializes the flow to the \c flowMap. The \c flowMap should
+ /// contain a feasible flow, ie. in each node excluding the source
+ /// and the target the incoming flow should be equal to the
+ /// outgoing flow.
+ /// \return %False when the given flowMap does not contain
+ /// feasible flow.
+ template <typename FlowMap>
+ bool checkedFlowInit(const FlowMap& flowMap) {
+ createStructures();
+
+ _dt->clear();
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ _dt->makeTree(n);
+ _dt->addCost(n, _max_value);
+ }
+
+ for (EdgeIt e(_graph); e != INVALID; ++e) {
+ _flow->set(e, flowMap[e]);
+ }
+ for (NodeIt it(_graph); it != INVALID; ++it) {
+ if (it == _source || it == _target) continue;
+ Value outFlow = 0;
+ for (OutEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
+ outFlow += (*_flow)[jt];
+ }
+ Value inFlow = 0;
+ for (InEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
+ inFlow += (*_flow)[jt];
+ }
+ if (_tolerance.different(outFlow, inFlow)) {
+ return false;
+ }
+ }
+ for (EdgeIt it(_graph); it != INVALID; ++it) {
+ if (_tolerance.less((*_flow)[it], 0)) return false;
+ if (_tolerance.less((*_capacity)[it], (*_flow)[it])) return false;
+ }
+ _flow_value = 0;
+ for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
+ _flow_value += (*_flow)[jt];
+ }
+ for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
+ _flow_value -= (*_flow)[jt];
+ }
+ return true;
+ }
+
+ /// \brief Executes the algorithm
+ ///
+ /// It runs augmenting phases by adding blocking flow until the
+ /// optimal solution is reached.
+ void start() {
+ while (augment());
+ }
+
+ /// \brief Augments the flow with a blocking flow on a layered
+ /// graph.
+ ///
+ /// This function builds a layered graph and then find a blocking
+ /// flow on this graph. The number of the levels in the layered
+ /// graph is strictly increasing in each augmenting phase
+ /// therefore the number of the augmentings is at most \f$ n-1
+ /// \f$. The length of each phase is at most \f$ O(m \log(n))
+ /// \f$, that the overall time complexity is \f$ O(nm \log(n)) \f$.
+ /// \return %False when there is not residual path between the
+ /// source and the target so the current flow is a feasible and
+ /// optimal solution.
+ bool augment() {
+ Node n;
+
+ if (createLayeredGraph()) {
+
+ Timer bf_timer;
+ initEdges();
+
+ n = _dt->findRoot(_source);
+ while (true) {
+ Edge e;
+ if (n == _target) {
+ augmentPath();
+ } else if (!advance(n)) {
+ if (n != _source) {
+ retreat(n);
+ } else {
+ break;
+ }
+ }
+ n = _dt->findRoot(_source);
+ }
+ extractTrees();
+
+ return true;
+ } else {
+ return false;
+ }
+ }
+
+ /// \brief runs the algorithm.
+ ///
+ /// It is just a shorthand for:
+ ///
+ ///\code
+ /// ek.init();
+ /// ek.start();
+ ///\endcode
+ void run() {
+ init();
+ start();
+ }
+
+ /// @}
+
+ /// \name Query Functions
+ /// The result of the %Dijkstra algorithm can be obtained using these
+ /// functions.\n
+ /// Before the use of these functions,
+ /// either run() or start() must be called.
+
+ ///@{
+
+ /// \brief Returns the value of the maximum flow.
+ ///
+ /// Returns the value of the maximum flow by returning the excess
+ /// of the target node \c t. This value equals to the value of
+ /// the maximum flow already after the first phase.
+ Value flowValue() const {
+ return _flow_value;
+ }
+
+
+ /// \brief Returns the flow on the edge.
+ ///
+ /// Sets the \c flowMap to the flow on the edges. This method can
+ /// be called after the second phase of algorithm.
+ Value flow(const Edge& edge) const {
+ return (*_flow)[edge];
+ }
+
+ /// \brief Returns true when the node is on the source side of minimum cut.
+ ///
+
+ /// Returns true when the node is on the source side of minimum
+ /// cut. This method can be called both after running \ref
+ /// startFirstPhase() and \ref startSecondPhase().
+ bool minCut(const Node& node) const {
+ return (*_level)[node] == -2;
+ }
+
+ /// \brief Returns a minimum value cut.
+ ///
+ /// Sets \c cut to the characteristic vector of a minimum value cut
+ /// It simply calls the minMinCut member.
+ /// \retval cut Write node bool map.
+ template <typename CutMap>
+ void minCutMap(CutMap& cutMap) const {
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ cutMap.set(n, (*_level)[n] == -2);
+ }
+ cutMap.set(_source, true);
+ }
+
+ /// @}
+
+ };
+}
+
+#endif
Added: lemon/trunk/lemon/dynamic_tree.h
==============================================================================
--- (empty file)
+++ lemon/trunk/lemon/dynamic_tree.h Sat Nov 17 21:58:11 2007
@@ -0,0 +1,1076 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+#ifndef LEMON_DYNAMIC_TREE_H
+#define LEMON_DYNAMIC_TREE_H
+
+/// \ingroup auxdata
+/// \file
+/// \brief The dynamic tree data structure of Sleator and Tarjan.
+///
+
+#include <vector>
+#include <limits>
+#include <lemon/tolerance.h>
+
+namespace lemon {
+
+ /// \ingroup auxdata
+ ///
+ /// \brief The dynamic tree data structure of Sleator and Tarjan.
+ ///
+ /// This class provides an implementation of the dynamic tree data
+ /// structure for maintaining a set of node-disjoint rooted
+ /// trees. Each item has an associated value, and the item with
+ /// minimum value can be find in \f$O(\log(n)\f$ on the path from a
+ /// node to the its root, and the items on such path can be
+ /// increased or decreased globally with a certain value in the same
+ /// running time. We regard a tree edge as directed toward the root,
+ /// that is from child to parent. Its structure can be modified by
+ /// two basic operations: \e link(v,w) adds an edge between a root v
+ /// and a node w in a different component; \e cut(v) removes the
+ /// edge between v and its parent.
+ ///
+ /// \param _Value The value type of the items.
+ /// \param _ItemIntMap Converts item type of node to integer.
+ /// \param _Tolerance The tolerance class to handle computation
+ /// problems.
+ /// \param _enableSize If true then the data structre manatain the
+ /// size of each tree. The feature is used in \ref GoldbergTarjan
+ /// algorithm. The default value is true.
+ ///
+ /// \author Hamori Tamas
+#ifdef DOXYGEN
+ template <typename _Value, typename _ItemIntMap,
+ typename _Tolerance, bool _enableSize>
+#else
+ template <typename _Value, typename _ItemIntMap,
+ typename _Tolerance = lemon::Tolerance<_Value>,
+ bool _enableSize = true>
+#endif
+ class DynamicTree {
+ public:
+ /// \brief The integer map on the items.
+ typedef _ItemIntMap ItemIntMap;
+ /// \brief The item type of nodes.
+ typedef typename ItemIntMap::Key Item;
+ /// \brief The value type of the algorithms.
+ typedef _Value Value;
+ /// \brief The tolerance used by the algorithm.
+ typedef _Tolerance Tolerance;
+
+ private:
+
+ class ItemData;
+
+ std::vector<ItemData> _data;
+ ItemIntMap &_iim;
+ Value _max_value;
+ Tolerance _tolerance;
+
+ public:
+ /// \brief The constructor of the class.
+ ///
+ /// \param iim The integer map on the items.
+ /// \param tolerance Tolerance class.
+ DynamicTree(ItemIntMap &iim, const Tolerance& tolerance = Tolerance())
+ : _iim(iim), _max_value(std::numeric_limits<Value>::max()),
+ _tolerance(tolerance) {}
+
+ ~DynamicTree() {}
+
+ /// \brief Clears the data structure
+ ///
+ /// Clears the data structure
+ void clear() {
+ _data.clear();
+ }
+
+ /// \brief Sets the tolerance used by algorithm.
+ ///
+ /// Sets the tolerance used by algorithm.
+ void tolerance(const Tolerance& tolerance) const {
+ _tolerance = tolerance;
+ return *this;
+ }
+
+ /// \brief Returns the tolerance used by algorithm.
+ ///
+ /// Returns the tolerance used by algorithm.
+ const Tolerance& tolerance() const {
+ return tolerance;
+ }
+
+ /// \brief Create a new tree containing a single node with cost zero.
+ void makeTree(const Item &item) {
+ _data[makePath(item)].successor = -1;
+ }
+
+ /// \brief Return the root of the tree containing node with itemtype
+ /// \e item.
+ Item findRoot(const Item &item) {
+ return _data[findTail(expose(_iim[item]))].id;
+ }
+
+ /// \brief Return the the value of nodes in the tree containing
+ /// node with itemtype \e item.
+ int findSize(const Item &item) {
+ return _data[expose(_iim[item])].size;
+ }
+
+ /// \brief Return the minimum cost containing node.
+ ///
+ /// Return into \e d the minimum cost on the tree path from \e item
+ /// to findRoot(item). Return the last item (closest to its root)
+ /// on this path of the minimum cost.
+ Item findCost(const Item &item, Value& d){
+ return _data[findPathCost(expose(_iim[item]),d)].id;
+ }
+
+ /// \brief Add \e x value to the cost of every node on the path from
+ /// \e item to findRoot(item).
+ void addCost(const Item &item, Value x){
+ addPathCost(expose(_iim[item]), x);
+ }
+
+ /// \brief Combine the trees containing nodes \e item1 and \e item2
+ /// by adding an edge from \e item1 \e item2.
+ ///
+ /// This operation assumes that \e item1 is root and \e item2 is in
+ /// a different tree.
+ void link(const Item &item1, const Item &item2){
+ int v = _iim[item1];
+ int w = _iim[item2];
+ int p = expose(w);
+ join(-1, expose(v), p);
+ _data[v].successor = -1;
+ _data[v].size += _data[p].size;
+
+ }
+
+ /// \brief Divide the tree containing node \e item into two trees by
+ /// deleting the edge out of \e item.
+ ///
+ /// This operation assumes that \e item is not a tree root.
+ void cut(const Item &item) {
+ int v = _iim[item];
+ int p, q;
+ expose(v);
+ split(p, v, q);
+ if (p != -1) {
+ _data[p].successor = v;
+ }
+ _data[v].size -= _data[q].size;
+ if (q != -1) {
+ _data[q].successor = _data[v].successor;
+ }
+ _data[v].successor = -1;
+ }
+
+ ///\brief
+ Item parent(const Item &item){
+ return _data[_iim[item].p].id;
+ }
+
+ ///\brief Return the upper bound of the costs.
+ Value maxValue() const {
+ return _max_value;
+ }
+
+ private:
+
+ int makePath(const Item &item) {
+ _iim.set(item, _data.size());
+ ItemData v(item);
+ _data.push_back(v);
+ return _iim[item];
+ }
+
+ int findPath(int v){
+ splay(v);
+ return v;
+ }
+
+ int findPathCost(int p, Value &d){
+ while ((_data[p].right != -1 &&
+ !_tolerance.less(0, _data[_data[p].right].dmin)) ||
+ (_data[p].left != -1 && _tolerance.less(0, _data[p].dcost))) {
+ if (_data[p].right != -1 &&
+ !_tolerance.less(0, _data[_data[p].right].dmin)) {
+ p = _data[p].right;
+ } else if (_data[p].left != -1 &&
+ !_tolerance.less(0, _data[_data[p].left].dmin)){
+ p = _data[p].left;
+ }
+ }
+ splay(p);
+ d = _data[p].dmin;
+ return p;
+ }
+
+ int findTail(int p){
+ while (_data[p].right != -1) {
+ p = _data[p].right;
+ }
+ splay(p);
+ return p;
+ }
+
+ void addPathCost(int p, Value x){
+ if (!_tolerance.less(x, _max_value)) {
+ _data[p].dmin = x;_data[p].dcost = x;
+ } else if (!_tolerance.less(-x, _max_value)) {
+ _data[p].dmin = 0;
+ _data[p].dcost = 0;
+ } else {
+ _data[p].dmin += x;
+ }
+ }
+
+ void join(int p, int v, int q) {
+ Value min = _max_value;
+ Value pmin = _max_value;
+ Value vmin = _data[v].dmin;
+ Value qmin = _max_value;
+ if (p != -1){
+ pmin = _data[p].dmin;
+ }
+ if (q != -1){
+ qmin = _data[q].dmin;
+ }
+
+ if (_tolerance.less(vmin, qmin)) {
+ if (_tolerance.less(vmin,pmin)) {
+ min = vmin;
+ } else {
+ min = pmin;
+ }
+ } else if (_tolerance.less(qmin,pmin)) {
+ min = qmin;
+ } else {
+ min = pmin;
+ }
+
+ if (p != -1){
+ _data[p].parent = v;
+ _data[p].dmin -= min;
+ }
+ if (q!=-1){
+ _data[q].parent = v;
+ if (_tolerance.less(_data[q].dmin,_max_value)) {
+ _data[q].dmin -= min;
+ }
+ }
+ _data[v].left = p;
+ _data[v].right = q;
+ if (_tolerance.less(min,_max_value)) {
+ _data[v].dcost = _data[v].dmin - min;
+ }
+ _data[v].dmin = min;
+ }
+
+ void split(int &p, int v, int &q){
+ splay(v);
+ p = -1;
+ if (_data[v].left != -1){
+ p = _data[v].left;
+ _data[p].dmin += _data[v].dmin;
+ _data[p].parent = -1;
+ _data[v].left = -1;
+ }
+ q = -1;
+ if (_data[v].right != -1) {
+ q=_data[v].right;
+ if (_tolerance.less(_data[q].dmin, _max_value)) {
+ _data[q].dmin += _data[v].dmin;
+ }
+ _data[q].parent = -1;
+ _data[v].right = -1;
+ }
+ if (_tolerance.less(_data[v].dcost, _max_value)) {
+ _data[v].dmin += _data[v].dcost;
+ _data[v].dcost = 0;
+ } else {
+ _data[v].dmin = _data[v].dcost;
+ }
+ }
+
+ int expose(int v) {
+ int p, q, r, w;
+ p = -1;
+ while (v != -1) {
+ w = _data[findPath(v)].successor;
+ split(q, v, r);
+ if (q != -1) {
+ _data[q].successor = v;
+ }
+ join(p, v, r);
+ p = v;
+ v = w;
+ }
+ _data[p].successor = -1;
+ return p;
+ }
+
+ void splay(int v) {
+ while (_data[v].parent != -1) {
+ if (v == _data[_data[v].parent].left) {
+ if (_data[_data[v].parent].parent == -1) {
+ zig(v);
+ } else {
+ if (_data[v].parent == _data[_data[_data[v].parent].parent].left) {
+ zig(_data[v].parent);
+ zig(v);
+ } else {
+ zig(v);
+ zag(v);
+ }
+ }
+ } else {
+ if (_data[_data[v].parent].parent == -1) {
+ zag(v);
+ } else {
+ if (_data[v].parent == _data[_data[_data[v].parent].parent].left) {
+ zag(v);
+ zig(v);
+ } else {
+ zag(_data[v].parent);
+ zag(v);
+ }
+ }
+ }
+ }
+ }
+
+
+ void zig(int v) {
+ Value min = _data[_data[v].parent].dmin;
+ int a = _data[v].parent;
+
+ Value aa = _data[a].dcost;
+ if (_tolerance.less(aa, _max_value)) {
+ aa+= min;
+ }
+
+
+ int b = v;
+ Value ab = min + _data[b].dmin;
+ Value bb = _data[b].dcost;
+ if (_tolerance.less(bb, _max_value)) {
+ bb+= ab;
+ }
+
+ int c = -1;
+ Value cc = _max_value;
+ if (_data[a].right != -1) {
+ c = _data[a].right;
+ cc = _data[c].dmin;
+ if (_tolerance.less(cc, _max_value)) {
+ cc+=min;
+ }
+ }
+
+ int d = -1;
+ Value dd = _max_value;
+ if (_data[v].left != -1){
+ d = _data[v].left;
+ dd = ab + _data[d].dmin;
+ }
+
+ int e = -1;
+ Value ee = _max_value;
+ if (_data[v].right != -1) {
+ e = _data[v].right;
+ ee = ab + _data[e].dmin;
+ }
+
+ Value min2;
+ if (_tolerance.less(0, _data[b].dmin) ||
+ (e != -1 && !_tolerance.less(0, _data[e].dmin))) {
+ min2 = min;
+ } else {
+ if (_tolerance.less(aa, cc)) {
+ if (_tolerance.less(aa, ee)) {
+ min2 = aa;
+ } else {
+ min2 = ee;
+ }
+ } else if (_tolerance.less(cc, ee)) {
+ min2 = cc;
+ } else {
+ min2 = ee;
+ }
+ }
+
+ _data[a].dcost = aa;
+ if (_tolerance.less(aa, _max_value)) {
+ _data[a].dcost -= min2;
+ }
+ _data[a].dmin = min2;
+ if (_tolerance.less(min2,_max_value)) {
+ _data[a].dmin -= min;
+ }
+ _data[a].size -= _data[b].size;
+ _data[b].dcost = bb;
+ if (_tolerance.less(bb, _max_value)) {
+ _data[b].dcost -= min;
+ }
+ _data[b].dmin = min;
+ _data[b].size += _data[a].size;
+ if (c != -1) {
+ _data[c].dmin = cc;
+ if (_tolerance.less(cc, _max_value)) {
+ _data[c].dmin -= min2;
+ }
+ }
+ if (d != -1) {
+ _data[d].dmin = dd - min;
+ _data[a].size += _data[d].size;
+ _data[b].size -= _data[d].size;
+ }
+ if (e != -1) {
+ _data[e].dmin = ee - min2;
+ }
+
+ int w = _data[v].parent;
+ _data[v].successor = _data[w].successor;
+ _data[w].successor = -1;
+ _data[v].parent = _data[w].parent;
+ _data[w].parent = v;
+ _data[w].left = _data[v].right;
+ _data[v].right = w;
+ if (_data[v].parent != -1){
+ if (_data[_data[v].parent].right == w) {
+ _data[_data[v].parent].right = v;
+ } else {
+ _data[_data[v].parent].left = v;
+ }
+ }
+ if (_data[w].left != -1){
+ _data[_data[w].left].parent = w;
+ }
+ }
+
+
+ void zag(int v) {
+
+ Value min = _data[_data[v].parent].dmin;
+
+ int a = _data[v].parent;
+ Value aa = _data[a].dcost;
+ if (_tolerance.less(aa, _max_value)) {
+ aa += min;
+ }
+
+ int b = v;
+ Value ab = min + _data[b].dmin;
+ Value bb = _data[b].dcost;
+ if (_tolerance.less(bb, _max_value)) {
+ bb += ab;
+ }
+
+ int c = -1;
+ Value cc = _max_value;
+ if (_data[a].left != -1){
+ c = _data[a].left;
+ cc = min + _data[c].dmin;
+ }
+
+ int d = -1;
+ Value dd = _max_value;
+ if (_data[v].right!=-1) {
+ d = _data[v].right;
+ dd = _data[d].dmin;
+ if (_tolerance.less(dd, _max_value)) {
+ dd += ab;
+ }
+ }
+
+ int e = -1;
+ Value ee = _max_value;
+ if (_data[v].left != -1){
+ e = _data[v].left;
+ ee = ab + _data[e].dmin;
+ }
+
+ Value min2;
+ if (_tolerance.less(0, _data[b].dmin) ||
+ (e != -1 && !_tolerance.less(0, _data[e].dmin))) {
+ min2 = min;
+ } else {
+ if (_tolerance.less(aa, cc)) {
+ if (_tolerance.less(aa, ee)) {
+ min2 = aa;
+ } else {
+ min2 = ee;
+ }
+ } else if (_tolerance.less(cc, ee)) {
+ min2 = cc;
+ } else {
+ min2 = ee;
+ }
+ }
+ _data[a].dcost = aa;
+ if (_tolerance.less(aa, _max_value)) {
+ _data[a].dcost -= min2;
+ }
+ _data[a].dmin = min2;
+ if (_tolerance.less(min2, _max_value)) {
+ _data[a].dmin -= min;
+ }
+ _data[a].size -= _data[b].size;
+ _data[b].dcost = bb;
+ if (_tolerance.less(bb, _max_value)) {
+ _data[b].dcost -= min;
+ }
+ _data[b].dmin = min;
+ _data[b].size += _data[a].size;
+ if (c != -1) {
+ _data[c].dmin = cc - min2;
+ }
+ if (d != -1) {
+ _data[d].dmin = dd;
+ _data[a].size += _data[d].size;
+ _data[b].size -= _data[d].size;
+ if (_tolerance.less(dd, _max_value)) {
+ _data[d].dmin -= min;
+ }
+ }
+ if (e != -1) {
+ _data[e].dmin = ee - min2;
+ }
+
+ int w = _data[v].parent;
+ _data[v].successor = _data[w].successor;
+ _data[w].successor = -1;
+ _data[v].parent = _data[w].parent;
+ _data[w].parent = v;
+ _data[w].right = _data[v].left;
+ _data[v].left = w;
+ if (_data[v].parent != -1){
+ if (_data[_data[v].parent].left == w) {
+ _data[_data[v].parent].left = v;
+ } else {
+ _data[_data[v].parent].right = v;
+ }
+ }
+ if (_data[w].right != -1){
+ _data[_data[w].right].parent = w;
+ }
+ }
+
+ private:
+
+ class ItemData {
+ public:
+ Item id;
+ int size;
+ int successor;
+ int parent;
+ int left;
+ int right;
+ Value dmin;
+ Value dcost;
+
+ public:
+ ItemData(const Item &item)
+ : id(item), size(1), successor(), parent(-1),
+ left(-1), right(-1), dmin(0), dcost(0) {}
+ };
+
+ };
+
+ template <typename _Value, typename _ItemIntMap, typename _Tolerance>
+ class DynamicTree<_Value, _ItemIntMap, _Tolerance, false> {
+ public:
+ typedef _ItemIntMap ItemIntMap;
+ typedef typename ItemIntMap::Key Item;
+ typedef _Value Value;
+ typedef _Tolerance Tolerance;
+
+ private:
+
+ class ItemData;
+
+ std::vector<ItemData> _data;
+ ItemIntMap &_iim;
+ Value _max_value;
+ Tolerance _tolerance;
+
+ public:
+ DynamicTree(ItemIntMap &iim, const Tolerance& tolerance = Tolerance())
+ : _iim(iim), _max_value(std::numeric_limits<Value>::max()),
+ _tolerance(tolerance) {}
+
+ ~DynamicTree() {}
+
+ void clear() {
+ _data.clear();
+ }
+
+ void tolerance(const Tolerance& tolerance) const {
+ _tolerance = tolerance;
+ return *this;
+ }
+
+ const Tolerance& tolerance() const {
+ return tolerance;
+ }
+
+ void makeTree(const Item &item) {
+ _data[makePath(item)].successor = -1;
+ }
+
+ Item findRoot(const Item &item) {
+ return _data[findTail(expose(_iim[item]))].id;
+ }
+
+ Item findCost(const Item &item, Value& d){
+ return _data[findPathCost(expose(_iim[item]),d)].id;
+ }
+
+ void addCost(const Item &item, Value x){
+ addPathCost(expose(_iim[item]), x);
+ }
+
+ void link(const Item &item1, const Item &item2){
+ int v = _iim[item1];
+ int w = _iim[item2];
+ int p = expose(w);
+ join(-1, expose(v), p);
+ _data[v].successor = -1;
+ }
+
+ void cut(const Item &item) {
+ int v = _iim[item];
+ int p, q;
+ expose(v);
+ split(p, v, q);
+ if (p != -1) {
+ _data[p].successor = v;
+ }
+ if (q != -1) {
+ _data[q].successor = _data[v].successor;
+ }
+ _data[v].successor = -1;
+ }
+
+ Item parent(const Item &item){
+ return _data[_iim[item].p].id;
+ }
+
+ Value maxValue() const {
+ return _max_value;
+ }
+
+ private:
+
+ int makePath(const Item &item) {
+ _iim.set(item, _data.size());
+ ItemData v(item);
+ _data.push_back(v);
+ return _iim[item];
+ }
+
+ int findPath(int v){
+ splay(v);
+ return v;
+ }
+
+ int findPathCost(int p, Value &d){
+ while ((_data[p].right != -1 &&
+ !_tolerance.less(0, _data[_data[p].right].dmin)) ||
+ (_data[p].left != -1 && _tolerance.less(0, _data[p].dcost))) {
+ if (_data[p].right != -1 &&
+ !_tolerance.less(0, _data[_data[p].right].dmin)) {
+ p = _data[p].right;
+ } else if (_data[p].left != -1 &&
+ !_tolerance.less(0, _data[_data[p].left].dmin)){
+ p = _data[p].left;
+ }
+ }
+ splay(p);
+ d = _data[p].dmin;
+ return p;
+ }
+
+ int findTail(int p){
+ while (_data[p].right != -1) {
+ p = _data[p].right;
+ }
+ splay(p);
+ return p;
+ }
+
+ void addPathCost(int p, Value x){
+ if (!_tolerance.less(x, _max_value)) {
+ _data[p].dmin = x;_data[p].dcost = x;
+ } else if (!_tolerance.less(-x, _max_value)) {
+ _data[p].dmin = 0;
+ _data[p].dcost = 0;
+ } else {
+ _data[p].dmin += x;
+ }
+ }
+
+ void join(int p, int v, int q) {
+ Value min = _max_value;
+ Value pmin = _max_value;
+ Value vmin = _data[v].dmin;
+ Value qmin = _max_value;
+ if (p != -1){
+ pmin = _data[p].dmin;
+ }
+ if (q != -1){
+ qmin = _data[q].dmin;
+ }
+
+ if (_tolerance.less(vmin, qmin)) {
+ if (_tolerance.less(vmin,pmin)) {
+ min = vmin;
+ } else {
+ min = pmin;
+ }
+ } else if (_tolerance.less(qmin,pmin)) {
+ min = qmin;
+ } else {
+ min = pmin;
+ }
+
+ if (p != -1){
+ _data[p].parent = v;
+ _data[p].dmin -= min;
+ }
+ if (q!=-1){
+ _data[q].parent = v;
+ if (_tolerance.less(_data[q].dmin,_max_value)) {
+ _data[q].dmin -= min;
+ }
+ }
+ _data[v].left = p;
+ _data[v].right = q;
+ if (_tolerance.less(min,_max_value)) {
+ _data[v].dcost = _data[v].dmin - min;
+ }
+ _data[v].dmin = min;
+ }
+
+ void split(int &p, int v, int &q){
+ splay(v);
+ p = -1;
+ if (_data[v].left != -1){
+ p = _data[v].left;
+ _data[p].dmin += _data[v].dmin;
+ _data[p].parent = -1;
+ _data[v].left = -1;
+ }
+ q = -1;
+ if (_data[v].right != -1) {
+ q=_data[v].right;
+ if (_tolerance.less(_data[q].dmin, _max_value)) {
+ _data[q].dmin += _data[v].dmin;
+ }
+ _data[q].parent = -1;
+ _data[v].right = -1;
+ }
+ if (_tolerance.less(_data[v].dcost, _max_value)) {
+ _data[v].dmin += _data[v].dcost;
+ _data[v].dcost = 0;
+ } else {
+ _data[v].dmin = _data[v].dcost;
+ }
+ }
+
+ int expose(int v) {
+ int p, q, r, w;
+ p = -1;
+ while (v != -1) {
+ w = _data[findPath(v)].successor;
+ split(q, v, r);
+ if (q != -1) {
+ _data[q].successor = v;
+ }
+ join(p, v, r);
+ p = v;
+ v = w;
+ }
+ _data[p].successor = -1;
+ return p;
+ }
+
+ void splay(int v) {
+ while (_data[v].parent != -1) {
+ if (v == _data[_data[v].parent].left) {
+ if (_data[_data[v].parent].parent == -1) {
+ zig(v);
+ } else {
+ if (_data[v].parent == _data[_data[_data[v].parent].parent].left) {
+ zig(_data[v].parent);
+ zig(v);
+ } else {
+ zig(v);
+ zag(v);
+ }
+ }
+ } else {
+ if (_data[_data[v].parent].parent == -1) {
+ zag(v);
+ } else {
+ if (_data[v].parent == _data[_data[_data[v].parent].parent].left) {
+ zag(v);
+ zig(v);
+ } else {
+ zag(_data[v].parent);
+ zag(v);
+ }
+ }
+ }
+ }
+ }
+
+
+ void zig(int v) {
+ Value min = _data[_data[v].parent].dmin;
+ int a = _data[v].parent;
+
+ Value aa = _data[a].dcost;
+ if (_tolerance.less(aa, _max_value)) {
+ aa+= min;
+ }
+
+
+ int b = v;
+ Value ab = min + _data[b].dmin;
+ Value bb = _data[b].dcost;
+ if (_tolerance.less(bb, _max_value)) {
+ bb+= ab;
+ }
+
+ int c = -1;
+ Value cc = _max_value;
+ if (_data[a].right != -1) {
+ c = _data[a].right;
+ cc = _data[c].dmin;
+ if (_tolerance.less(cc, _max_value)) {
+ cc+=min;
+ }
+ }
+
+ int d = -1;
+ Value dd = _max_value;
+ if (_data[v].left != -1){
+ d = _data[v].left;
+ dd = ab + _data[d].dmin;
+ }
+
+ int e = -1;
+ Value ee = _max_value;
+ if (_data[v].right != -1) {
+ e = _data[v].right;
+ ee = ab + _data[e].dmin;
+ }
+
+ Value min2;
+ if (_tolerance.less(0, _data[b].dmin) ||
+ (e != -1 && !_tolerance.less(0, _data[e].dmin))) {
+ min2 = min;
+ } else {
+ if (_tolerance.less(aa, cc)) {
+ if (_tolerance.less(aa, ee)) {
+ min2 = aa;
+ } else {
+ min2 = ee;
+ }
+ } else if (_tolerance.less(cc, ee)) {
+ min2 = cc;
+ } else {
+ min2 = ee;
+ }
+ }
+
+ _data[a].dcost = aa;
+ if (_tolerance.less(aa, _max_value)) {
+ _data[a].dcost -= min2;
+ }
+ _data[a].dmin = min2;
+ if (_tolerance.less(min2,_max_value)) {
+ _data[a].dmin -= min;
+ }
+ _data[b].dcost = bb;
+ if (_tolerance.less(bb, _max_value)) {
+ _data[b].dcost -= min;
+ }
+ _data[b].dmin = min;
+ if (c != -1) {
+ _data[c].dmin = cc;
+ if (_tolerance.less(cc, _max_value)) {
+ _data[c].dmin -= min2;
+ }
+ }
+ if (d != -1) {
+ _data[d].dmin = dd - min;
+ }
+ if (e != -1) {
+ _data[e].dmin = ee - min2;
+ }
+
+ int w = _data[v].parent;
+ _data[v].successor = _data[w].successor;
+ _data[w].successor = -1;
+ _data[v].parent = _data[w].parent;
+ _data[w].parent = v;
+ _data[w].left = _data[v].right;
+ _data[v].right = w;
+ if (_data[v].parent != -1){
+ if (_data[_data[v].parent].right == w) {
+ _data[_data[v].parent].right = v;
+ } else {
+ _data[_data[v].parent].left = v;
+ }
+ }
+ if (_data[w].left != -1){
+ _data[_data[w].left].parent = w;
+ }
+ }
+
+
+ void zag(int v) {
+
+ Value min = _data[_data[v].parent].dmin;
+
+ int a = _data[v].parent;
+ Value aa = _data[a].dcost;
+ if (_tolerance.less(aa, _max_value)) {
+ aa += min;
+ }
+
+ int b = v;
+ Value ab = min + _data[b].dmin;
+ Value bb = _data[b].dcost;
+ if (_tolerance.less(bb, _max_value)) {
+ bb += ab;
+ }
+
+ int c = -1;
+ Value cc = _max_value;
+ if (_data[a].left != -1){
+ c = _data[a].left;
+ cc = min + _data[c].dmin;
+ }
+
+ int d = -1;
+ Value dd = _max_value;
+ if (_data[v].right!=-1) {
+ d = _data[v].right;
+ dd = _data[d].dmin;
+ if (_tolerance.less(dd, _max_value)) {
+ dd += ab;
+ }
+ }
+
+ int e = -1;
+ Value ee = _max_value;
+ if (_data[v].left != -1){
+ e = _data[v].left;
+ ee = ab + _data[e].dmin;
+ }
+
+ Value min2;
+ if (_tolerance.less(0, _data[b].dmin) ||
+ (e != -1 && !_tolerance.less(0, _data[e].dmin))) {
+ min2 = min;
+ } else {
+ if (_tolerance.less(aa, cc)) {
+ if (_tolerance.less(aa, ee)) {
+ min2 = aa;
+ } else {
+ min2 = ee;
+ }
+ } else if (_tolerance.less(cc, ee)) {
+ min2 = cc;
+ } else {
+ min2 = ee;
+ }
+ }
+ _data[a].dcost = aa;
+ if (_tolerance.less(aa, _max_value)) {
+ _data[a].dcost -= min2;
+ }
+ _data[a].dmin = min2;
+ if (_tolerance.less(min2, _max_value)) {
+ _data[a].dmin -= min;
+ }
+ _data[b].dcost = bb;
+ if (_tolerance.less(bb, _max_value)) {
+ _data[b].dcost -= min;
+ }
+ _data[b].dmin = min;
+ if (c != -1) {
+ _data[c].dmin = cc - min2;
+ }
+ if (d != -1) {
+ _data[d].dmin = dd;
+ if (_tolerance.less(dd, _max_value)) {
+ _data[d].dmin -= min;
+ }
+ }
+ if (e != -1) {
+ _data[e].dmin = ee - min2;
+ }
+
+ int w = _data[v].parent;
+ _data[v].successor = _data[w].successor;
+ _data[w].successor = -1;
+ _data[v].parent = _data[w].parent;
+ _data[w].parent = v;
+ _data[w].right = _data[v].left;
+ _data[v].left = w;
+ if (_data[v].parent != -1){
+ if (_data[_data[v].parent].left == w) {
+ _data[_data[v].parent].left = v;
+ } else {
+ _data[_data[v].parent].right = v;
+ }
+ }
+ if (_data[w].right != -1){
+ _data[_data[w].right].parent = w;
+ }
+ }
+
+ private:
+
+ class ItemData {
+ public:
+ Item id;
+ int successor;
+ int parent;
+ int left;
+ int right;
+ Value dmin;
+ Value dcost;
+
+ public:
+ ItemData(const Item &item)
+ : id(item), successor(), parent(-1),
+ left(-1), right(-1), dmin(0), dcost(0) {}
+ };
+
+ };
+
+}
+
+#endif
Modified: lemon/trunk/lemon/edmonds_karp.h
==============================================================================
--- lemon/trunk/lemon/edmonds_karp.h (original)
+++ lemon/trunk/lemon/edmonds_karp.h Sat Nov 17 21:58:11 2007
@@ -23,47 +23,95 @@
/// \ingroup max_flow
/// \brief Implementation of the Edmonds-Karp algorithm.
-#include <lemon/graph_adaptor.h>
#include <lemon/tolerance.h>
-#include <lemon/bfs.h>
+#include <vector>
namespace lemon {
+ /// \brief Default traits class of EdmondsKarp class.
+ ///
+ /// Default traits class of EdmondsKarp class.
+ /// \param _Graph Graph type.
+ /// \param _CapacityMap Type of capacity map.
+ template <typename _Graph, typename _CapacityMap>
+ struct EdmondsKarpDefaultTraits {
+
+ /// \brief The graph type the algorithm runs on.
+ typedef _Graph Graph;
+
+ /// \brief The type of the map that stores the edge capacities.
+ ///
+ /// The type of the map that stores the edge capacities.
+ /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
+ typedef _CapacityMap CapacityMap;
+
+ /// \brief The type of the length of the edges.
+ typedef typename CapacityMap::Value Value;
+
+ /// \brief The map type that stores the flow values.
+ ///
+ /// The map type that stores the flow values.
+ /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
+ typedef typename Graph::template EdgeMap<Value> FlowMap;
+
+ /// \brief Instantiates a FlowMap.
+ ///
+ /// This function instantiates a \ref FlowMap.
+ /// \param graph The graph, to which we would like to define the flow map.
+ static FlowMap* createFlowMap(const Graph& graph) {
+ return new FlowMap(graph);
+ }
+
+ /// \brief The tolerance used by the algorithm
+ ///
+ /// The tolerance used by the algorithm to handle inexact computation.
+ typedef Tolerance<Value> Tolerance;
+
+ };
+
/// \ingroup max_flow
+ ///
/// \brief Edmonds-Karp algorithms class.
///
/// This class provides an implementation of the \e Edmonds-Karp \e
/// algorithm producing a flow of maximum value in a directed
- /// graph. The Edmonds-Karp algorithm is slower than the Preflow algorithm
- /// but it has an advantage of the step-by-step execution control with
- /// feasible flow solutions. The \e source node, the \e target node, the \e
- /// capacity of the edges and the \e starting \e flow value of the
- /// edges should be passed to the algorithm through the
- /// constructor.
+ /// graphs. The Edmonds-Karp algorithm is slower than the Preflow
+ /// algorithm but it has an advantage of the step-by-step execution
+ /// control with feasible flow solutions. The \e source node, the \e
+ /// target node, the \e capacity of the edges and the \e starting \e
+ /// flow value of the edges should be passed to the algorithm
+ /// through the constructor.
///
- /// The time complexity of the algorithm is \f$ O(n * e^2) \f$ in
+ /// The time complexity of the algorithm is \f$ O(nm^2) \f$ in
/// worst case. Always try the preflow algorithm instead of this if
/// you just want to compute the optimal flow.
///
/// \param _Graph The directed graph type the algorithm runs on.
- /// \param _Number The number type of the capacities and the flow values.
/// \param _CapacityMap The capacity map type.
- /// \param _FlowMap The flow map type.
- /// \param _Tolerance The tolerance class to handle computation problems.
+ /// \param _Traits Traits class to set various data types used by
+ /// the algorithm. The default traits class is \ref
+ /// EdmondsKarpDefaultTraits. See \ref EdmondsKarpDefaultTraits for the
+ /// documentation of a Edmonds-Karp traits class.
///
/// \author Balazs Dezso
#ifdef DOXYGEN
- template <typename _Graph, typename _Number,
- typename _CapacityMap, typename _FlowMap, typename _Tolerance>
-#else
- template <typename _Graph, typename _Number,
- typename _CapacityMap = typename _Graph::template EdgeMap<_Number>,
- typename _FlowMap = typename _Graph::template EdgeMap<_Number>,
- typename _Tolerance = Tolerance<_Number> >
+ template <typename _Graph, typename _CapacityMap, typename _Traits>
+#else
+ template <typename _Graph,
+ typename _CapacityMap = typename _Graph::template EdgeMap<int>,
+ typename _Traits = EdmondsKarpDefaultTraits<_Graph, _CapacityMap> >
#endif
class EdmondsKarp {
public:
+ typedef _Traits Traits;
+ typedef typename Traits::Graph Graph;
+ typedef typename Traits::CapacityMap CapacityMap;
+ typedef typename Traits::Value Value;
+
+ typedef typename Traits::FlowMap FlowMap;
+ typedef typename Traits::Tolerance Tolerance;
+
/// \brief \ref Exception for the case when the source equals the target.
///
/// \ref Exception for the case when the source equals the target.
@@ -76,110 +124,239 @@
};
- /// \brief The graph type the algorithm runs on.
- typedef _Graph Graph;
- /// \brief The value type of the algorithms.
- typedef _Number Number;
- /// \brief The capacity map on the edges.
- typedef _CapacityMap CapacityMap;
- /// \brief The flow map on the edges.
- typedef _FlowMap FlowMap;
- /// \brief The tolerance used by the algorithm.
- typedef _Tolerance Tolerance;
+ private:
- typedef ResGraphAdaptor<Graph, Number, CapacityMap,
- FlowMap, Tolerance> ResGraph;
+ GRAPH_TYPEDEFS(typename Graph);
+ typedef typename Graph::template NodeMap<Edge> PredMap;
+
+ const Graph& _graph;
+ const CapacityMap* _capacity;
- private:
+ Node _source, _target;
+
+ FlowMap* _flow;
+ bool _local_flow;
- typedef typename Graph::Node Node;
- typedef typename Graph::Edge Edge;
+ PredMap* _pred;
+ std::vector<Node> _queue;
- typedef typename Graph::NodeIt NodeIt;
- typedef typename Graph::EdgeIt EdgeIt;
- typedef typename Graph::InEdgeIt InEdgeIt;
- typedef typename Graph::OutEdgeIt OutEdgeIt;
+ Tolerance _tolerance;
+ Value _flow_value;
+
+ void createStructures() {
+ if (!_flow) {
+ _flow = Traits::createFlowMap(_graph);
+ _local_flow = true;
+ }
+ if (!_pred) {
+ _pred = new PredMap(_graph);
+ }
+ _queue.resize(countNodes(_graph));
+ }
+
+ void destroyStructures() {
+ if (_local_flow) {
+ delete _flow;
+ }
+ if (_pred) {
+ delete _pred;
+ }
+ }
public:
+ ///\name Named template parameters
+
+ ///@{
+
+ template <typename _FlowMap>
+ struct DefFlowMapTraits : public Traits {
+ typedef _FlowMap FlowMap;
+ static FlowMap *createFlowMap(const Graph&) {
+ throw UninitializedParameter();
+ }
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting
+ /// FlowMap type
+ ///
+ /// \ref named-templ-param "Named parameter" for setting FlowMap
+ /// type
+ template <typename _FlowMap>
+ struct DefFlowMap
+ : public EdmondsKarp<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > {
+ typedef EdmondsKarp<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> >
+ Create;
+ };
+
+
+ /// @}
+
/// \brief The constructor of the class.
///
/// The constructor of the class.
/// \param graph The directed graph the algorithm runs on.
+ /// \param capacity The capacity of the edges.
/// \param source The source node.
/// \param target The target node.
- /// \param capacity The capacity of the edges.
- /// \param flow The flow of the edges.
- /// \param tolerance Tolerance class.
- EdmondsKarp(const Graph& graph, Node source, Node target,
- const CapacityMap& capacity, FlowMap& flow,
- const Tolerance& tolerance = Tolerance())
- : _graph(graph), _capacity(capacity), _flow(flow),
- _tolerance(tolerance), _resgraph(graph, capacity, flow, tolerance),
- _source(source), _target(target)
+ EdmondsKarp(const Graph& graph, const CapacityMap& capacity,
+ Node source, Node target)
+ : _graph(graph), _capacity(&capacity), _source(source), _target(target),
+ _flow(0), _local_flow(false), _pred(0), _tolerance(), _flow_value()
{
if (_source == _target) {
throw InvalidArgument();
}
}
+ /// \brief Destrcutor.
+ ///
+ /// Destructor.
+ ~EdmondsKarp() {
+ destroyStructures();
+ }
+
+ /// \brief Sets the capacity map.
+ ///
+ /// Sets the capacity map.
+ /// \return \c (*this)
+ EdmondsKarp& capacityMap(const CapacityMap& map) {
+ _capacity = ↦
+ return *this;
+ }
+
+ /// \brief Sets the flow map.
+ ///
+ /// Sets the flow map.
+ /// \return \c (*this)
+ EdmondsKarp& flowMap(FlowMap& map) {
+ if (_local_flow) {
+ delete _flow;
+ _local_flow = false;
+ }
+ _flow = ↦
+ return *this;
+ }
+
+ /// \brief Returns the flow map.
+ ///
+ /// \return The flow map.
+ const FlowMap& flowMap() {
+ return *_flow;
+ }
+
+ /// \brief Sets the source node.
+ ///
+ /// Sets the source node.
+ /// \return \c (*this)
+ EdmondsKarp& source(const Node& node) {
+ _source = node;
+ return *this;
+ }
+
+ /// \brief Sets the target node.
+ ///
+ /// Sets the target node.
+ /// \return \c (*this)
+ EdmondsKarp& target(const Node& node) {
+ _target = node;
+ return *this;
+ }
+
+ /// \brief Sets the tolerance used by algorithm.
+ ///
+ /// Sets the tolerance used by algorithm.
+ EdmondsKarp& tolerance(const Tolerance& tolerance) const {
+ _tolerance = tolerance;
+ return *this;
+ }
+
+ /// \brief Returns the tolerance used by algorithm.
+ ///
+ /// Returns the tolerance used by algorithm.
+ const Tolerance& tolerance() const {
+ return tolerance;
+ }
+
+ /// \name Execution control The simplest way to execute the
+ /// algorithm is to use the \c run() member functions.
+ /// \n
+ /// If you need more control on initial solution or
+ /// execution then you have to call one \ref init() function and then
+ /// the start() or multiple times the \c augment() member function.
+
+ ///@{
+
/// \brief Initializes the algorithm
///
/// It sets the flow to empty flow.
void init() {
+ createStructures();
for (EdgeIt it(_graph); it != INVALID; ++it) {
- _flow.set(it, 0);
+ _flow->set(it, 0);
}
- _value = 0;
+ _flow_value = 0;
}
/// \brief Initializes the algorithm
///
- /// If the flow map initially flow this let the flow map
- /// unchanged but the flow value will be set by the flow
- /// on the outedges from the source.
- void flowInit() {
- _value = 0;
+ /// Initializes the flow to the \c flowMap. The \c flowMap should
+ /// contain a feasible flow, ie. in each node excluding the source
+ /// and the target the incoming flow should be equal to the
+ /// outgoing flow.
+ template <typename FlowMap>
+ void flowInit(const FlowMap& flowMap) {
+ createStructures();
+ for (EdgeIt e(_graph); e != INVALID; ++e) {
+ _flow->set(e, flowMap[e]);
+ }
+ _flow_value = 0;
for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
- _value += _flow[jt];
+ _flow_value += (*_flow)[jt];
}
for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
- _value -= _flow[jt];
+ _flow_value -= (*_flow)[jt];
}
}
/// \brief Initializes the algorithm
///
- /// If the flow map initially flow this let the flow map
- /// unchanged but the flow value will be set by the flow
- /// on the outedges from the source. It also checks that
- /// the flow map really contains a flow.
- /// \return %True when the flow map really a flow.
- bool checkedFlowInit() {
- _value = 0;
- for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
- _value += _flow[jt];
- }
- for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
- _value -= _flow[jt];
+ /// Initializes the flow to the \c flowMap. The \c flowMap should
+ /// contain a feasible flow, ie. in each node excluding the source
+ /// and the target the incoming flow should be equal to the
+ /// outgoing flow.
+ /// \return %False when the given flowMap does not contain
+ /// feasible flow.
+ template <typename FlowMap>
+ bool checkedFlowInit(const FlowMap& flowMap) {
+ createStructures();
+ for (EdgeIt e(_graph); e != INVALID; ++e) {
+ _flow->set(e, flowMap[e]);
}
for (NodeIt it(_graph); it != INVALID; ++it) {
if (it == _source || it == _target) continue;
- Number outFlow = 0;
+ Value outFlow = 0;
for (OutEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
- outFlow += _flow[jt];
+ outFlow += (*_flow)[jt];
}
- Number inFlow = 0;
+ Value inFlow = 0;
for (InEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
- inFlow += _flow[jt];
+ inFlow += (*_flow)[jt];
}
if (_tolerance.different(outFlow, inFlow)) {
return false;
}
}
for (EdgeIt it(_graph); it != INVALID; ++it) {
- if (_tolerance.less(_flow[it], 0)) return false;
- if (_tolerance.less(_capacity[it], _flow[it])) return false;
+ if (_tolerance.less((*_flow)[it], 0)) return false;
+ if (_tolerance.less((*_capacity)[it], (*_flow)[it])) return false;
+ }
+ _flow_value = 0;
+ for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
+ _flow_value += (*_flow)[jt];
+ }
+ for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
+ _flow_value -= (*_flow)[jt];
}
return true;
}
@@ -195,32 +372,76 @@
/// \return %False when the augmenting is not success so the
/// current flow is a feasible and optimal solution.
bool augment() {
- typename Bfs<ResGraph>
- ::template DefDistMap<NullMap<Node, int> >
- ::Create bfs(_resgraph);
-
- NullMap<Node, int> distMap;
- bfs.distMap(distMap);
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ _pred->set(n, INVALID);
+ }
- bfs.init();
- bfs.addSource(_source);
- bfs.start(_target);
-
- if (!bfs.reached(_target)) {
- return false;
- }
- Number min = _resgraph.rescap(bfs.predEdge(_target));
- for (Node it = bfs.predNode(_target); it != _source;
- it = bfs.predNode(it)) {
- if (min > _resgraph.rescap(bfs.predEdge(it))) {
- min = _resgraph.rescap(bfs.predEdge(it));
- }
- }
- for (Node it = _target; it != _source; it = bfs.predNode(it)) {
- _resgraph.augment(bfs.predEdge(it), min);
+ int first = 0, last = 1;
+
+ _queue[0] = _source;
+ _pred->set(_source, OutEdgeIt(_graph, _source));
+
+ while (first != last && (*_pred)[_target] == INVALID) {
+ Node n = _queue[first++];
+
+ for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_capacity)[e] - (*_flow)[e];
+ Node t = _graph.target(e);
+ if (_tolerance.positive(rem) && (*_pred)[t] == INVALID) {
+ _pred->set(t, e);
+ _queue[last++] = t;
+ }
+ }
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_flow)[e];
+ Node t = _graph.source(e);
+ if (_tolerance.positive(rem) && (*_pred)[t] == INVALID) {
+ _pred->set(t, e);
+ _queue[last++] = t;
+ }
+ }
+ }
+
+ if ((*_pred)[_target] != INVALID) {
+ Node n = _target;
+ Edge e = (*_pred)[n];
+
+ Value prem = (*_capacity)[e] - (*_flow)[e];
+ n = _graph.source(e);
+ while (n != _source) {
+ e = (*_pred)[n];
+ if (_graph.target(e) == n) {
+ Value rem = (*_capacity)[e] - (*_flow)[e];
+ if (rem < prem) prem = rem;
+ n = _graph.source(e);
+ } else {
+ Value rem = (*_flow)[e];
+ if (rem < prem) prem = rem;
+ n = _graph.target(e);
+ }
+ }
+
+ n = _target;
+ e = (*_pred)[n];
+
+ _flow->set(e, (*_flow)[e] + prem);
+ n = _graph.source(e);
+ while (n != _source) {
+ e = (*_pred)[n];
+ if (_graph.target(e) == n) {
+ _flow->set(e, (*_flow)[e] + prem);
+ n = _graph.source(e);
+ } else {
+ _flow->set(e, (*_flow)[e] - prem);
+ n = _graph.target(e);
+ }
+ }
+
+ _flow_value += prem;
+ return true;
+ } else {
+ return false;
}
- _value += min;
- return true;
}
/// \brief Executes the algorithm
@@ -230,13 +451,6 @@
while (augment()) {}
}
- /// \brief Gives back the current flow value.
- ///
- /// Gives back the current flow _value.
- Number flowValue() const {
- return _value;
- }
-
/// \brief runs the algorithm.
///
/// It is just a shorthand for:
@@ -250,119 +464,59 @@
start();
}
- /// \brief Returns a minimum value cut.
- ///
- /// Sets \c cut to the characteristic vector of a minimum value cut
- /// It simply calls the minMinCut member.
- /// \retval cut Write node bool map.
- template <typename CutMap>
- void minCut(CutMap& cut) const {
- minMinCut(cut);
- }
+ /// @}
- /// \brief Returns the inclusionwise minimum of the minimum value cuts.
- ///
- /// Sets \c cut to the characteristic vector of the minimum value cut
- /// which is inclusionwise minimum. It is computed by processing a
- /// bfs from the source node \c source in the residual graph.
- /// \retval cut Write node bool map.
- template <typename CutMap>
- void minMinCut(CutMap& cut) const {
-
- typename Bfs<ResGraph>
- ::template DefDistMap<NullMap<Node, int> >
- ::template DefProcessedMap<CutMap>
- ::Create bfs(_resgraph);
-
- NullMap<Node, int> distMap;
- bfs.distMap(distMap);
-
- bfs.processedMap(cut);
-
- bfs.run(_source);
- }
+ /// \name Query Functions
+ /// The result of the %Dijkstra algorithm can be obtained using these
+ /// functions.\n
+ /// Before the use of these functions,
+ /// either run() or start() must be called.
+
+ ///@{
- /// \brief Returns the inclusionwise minimum of the minimum value cuts.
+ /// \brief Returns the value of the maximum flow.
///
- /// Sets \c cut to the characteristic vector of the minimum value cut
- /// which is inclusionwise minimum. It is computed by processing a
- /// bfs from the source node \c source in the residual graph.
- /// \retval cut Write node bool map.
- template <typename CutMap>
- void maxMinCut(CutMap& cut) const {
-
- typedef RevGraphAdaptor<const ResGraph> RevGraph;
-
- RevGraph revgraph(_resgraph);
-
- typename Bfs<RevGraph>
- ::template DefDistMap<NullMap<Node, int> >
- ::template DefPredMap<NullMap<Node, Edge> >
- ::template DefProcessedMap<NotWriteMap<CutMap> >
- ::Create bfs(revgraph);
-
- NullMap<Node, int> distMap;
- bfs.distMap(distMap);
-
- NullMap<Node, Edge> predMap;
- bfs.predMap(predMap);
-
- NotWriteMap<CutMap> notcut(cut);
- bfs.processedMap(notcut);
-
- bfs.run(_target);
+ /// Returns the value of the maximum flow by returning the excess
+ /// of the target node \c t. This value equals to the value of
+ /// the maximum flow already after the first phase.
+ Value flowValue() const {
+ return _flow_value;
}
- /// \brief Returns the source node.
- ///
- /// Returns the source node.
- ///
- Node source() const {
- return _source;
- }
- /// \brief Returns the target node.
+ /// \brief Returns the flow on the edge.
///
- /// Returns the target node.
- ///
- Node target() const {
- return _target;
+ /// Sets the \c flowMap to the flow on the edges. This method can
+ /// be called after the second phase of algorithm.
+ Value flow(const Edge& edge) const {
+ return (*_flow)[edge];
}
- /// \brief Returns a reference to capacity map.
+ /// \brief Returns true when the node is on the source side of minimum cut.
///
- /// Returns a reference to capacity map.
- ///
- const CapacityMap &capacityMap() const {
- return *_capacity;
- }
-
- /// \brief Returns a reference to flow map.
- ///
- /// Returns a reference to flow map.
- ///
- const FlowMap &flowMap() const {
- return *_flow;
+
+ /// Returns true when the node is on the source side of minimum
+ /// cut. This method can be called both after running \ref
+ /// startFirstPhase() and \ref startSecondPhase().
+ bool minCut(const Node& node) const {
+ return (*_pred)[node] != INVALID;
}
- /// \brief Returns the tolerance used by algorithm.
+ /// \brief Returns a minimum value cut.
///
- /// Returns the tolerance used by algorithm.
- const Tolerance& tolerance() const {
- return tolerance;
- }
-
- private:
-
- const Graph& _graph;
- const CapacityMap& _capacity;
- FlowMap& _flow;
- Tolerance _tolerance;
-
- ResGraph _resgraph;
- Node _source, _target;
- Number _value;
-
+ /// Sets \c cut to the characteristic vector of a minimum value cut
+ /// It simply calls the minMinCut member.
+ /// \retval cut Write node bool map.
+ template <typename CutMap>
+ void minCutMap(CutMap& cutMap) const {
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ cutMap.set(n, (*_pred)[n] != INVALID);
+ }
+ cutMap.set(_source, true);
+ }
+
+ /// @}
+
};
}
Added: lemon/trunk/lemon/goldberg_tarjan.h
==============================================================================
--- (empty file)
+++ lemon/trunk/lemon/goldberg_tarjan.h Sat Nov 17 21:58:11 2007
@@ -0,0 +1,1020 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_GOLDBERG_TARJAN_H
+#define LEMON_GOLDBERG_TARJAN_H
+
+#include <vector>
+#include <queue>
+
+#include <lemon/error.h>
+#include <lemon/bits/invalid.h>
+#include <lemon/tolerance.h>
+#include <lemon/maps.h>
+#include <lemon/graph_utils.h>
+#include <lemon/dynamic_tree.h>
+#include <limits>
+
+/// \file
+/// \ingroup max_flow
+/// \brief Implementation of the preflow algorithm.
+
+namespace lemon {
+
+ /// \brief Default traits class of GoldbergTarjan class.
+ ///
+ /// Default traits class of GoldbergTarjan class.
+ /// \param _Graph Graph type.
+ /// \param _CapacityMap Type of capacity map.
+ template <typename _Graph, typename _CapacityMap>
+ struct GoldbergTarjanDefaultTraits {
+
+ /// \brief The graph type the algorithm runs on.
+ typedef _Graph Graph;
+
+ /// \brief The type of the map that stores the edge capacities.
+ ///
+ /// The type of the map that stores the edge capacities.
+ /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
+ typedef _CapacityMap CapacityMap;
+
+ /// \brief The type of the length of the edges.
+ typedef typename CapacityMap::Value Value;
+
+ /// \brief The map type that stores the flow values.
+ ///
+ /// The map type that stores the flow values.
+ /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
+ typedef typename Graph::template EdgeMap<Value> FlowMap;
+
+ /// \brief Instantiates a FlowMap.
+ ///
+ /// This function instantiates a \ref FlowMap.
+ /// \param graph The graph, to which we would like to define the flow map.
+ static FlowMap* createFlowMap(const Graph& graph) {
+ return new FlowMap(graph);
+ }
+
+ /// \brief The eleavator type used by GoldbergTarjan algorithm.
+ ///
+ /// The elevator type used by GoldbergTarjan algorithm.
+ ///
+ /// \sa Elevator
+ /// \sa LinkedElevator
+ typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
+
+ /// \brief Instantiates an Elevator.
+ ///
+ /// This function instantiates a \ref Elevator.
+ /// \param graph The graph, to which we would like to define the elevator.
+ /// \param max_level The maximum level of the elevator.
+ static Elevator* createElevator(const Graph& graph, int max_level) {
+ return new Elevator(graph, max_level);
+ }
+
+ /// \brief The tolerance used by the algorithm
+ ///
+ /// The tolerance used by the algorithm to handle inexact computation.
+ typedef Tolerance<Value> Tolerance;
+
+ };
+
+ /// \ingroup max_flow
+ /// \brief Goldberg-Tarjan algorithms class.
+ ///
+ /// This class provides an implementation of the \e GoldbergTarjan
+ /// \e algorithm producing a flow of maximum value in a directed
+ /// graph. The GoldbergTarjan algorithm is a theoretical improvement
+ /// of the Goldberg's \ref Preflow "preflow" algorithm by using the \ref
+ /// DynamicTree "dynamic tree" data structure of Sleator and Tarjan.
+ ///
+ /// The original preflow algorithm with \e "highest label" or \e
+ /// FIFO heuristic has \f$O(n^3)\f$ time complexity. The current
+ /// algorithm improved this complexity to
+ /// \f$O(nm\log(\frac{n^2}{m}))\f$. The algorithm builds limited
+ /// size dynamic trees on the residual graph, which can be used to
+ /// make multilevel push operations and gives a better bound on the
+ /// number of non-saturating pushes.
+ ///
+ /// \param Graph The directed graph type the algorithm runs on.
+ /// \param CapacityMap The capacity map type.
+ /// \param _Traits Traits class to set various data types used by
+ /// the algorithm. The default traits class is \ref
+ /// GoldbergTarjanDefaultTraits. See \ref
+ /// GoldbergTarjanDefaultTraits for the documentation of a
+ /// Goldberg-Tarjan traits class.
+ ///
+ /// \author Tamas Hamori and Balazs Dezso
+#ifdef DOXYGEN
+ template <typename _Graph, typename _CapacityMap, typename _Traits>
+#else
+ template <typename _Graph,
+ typename _CapacityMap = typename _Graph::template EdgeMap<int>,
+ typename _Traits =
+ GoldbergTarjanDefaultTraits<_Graph, _CapacityMap> >
+#endif
+ class GoldbergTarjan {
+ public:
+
+ typedef _Traits Traits;
+ typedef typename Traits::Graph Graph;
+ typedef typename Traits::CapacityMap CapacityMap;
+ typedef typename Traits::Value Value;
+
+ typedef typename Traits::FlowMap FlowMap;
+ typedef typename Traits::Elevator Elevator;
+ typedef typename Traits::Tolerance Tolerance;
+
+ protected:
+
+ GRAPH_TYPEDEFS(typename Graph);
+
+ typedef typename Graph::template NodeMap<Node> NodeNodeMap;
+ typedef typename Graph::template NodeMap<int> IntNodeMap;
+
+ typedef typename Graph::template NodeMap<Edge> EdgeNodeMap;
+ typedef typename Graph::template EdgeMap<Edge> EdgeEdgeMap;
+
+ typedef typename std::vector<Node> VecNode;
+
+ typedef DynamicTree<Value,IntNodeMap,Tolerance> DynTree;
+
+ const Graph& _graph;
+ const CapacityMap* _capacity;
+ int _node_num; //the number of nodes of G
+
+ Node _source;
+ Node _target;
+
+ FlowMap* _flow;
+ bool _local_flow;
+
+ Elevator* _level;
+ bool _local_level;
+
+ typedef typename Graph::template NodeMap<Value> ExcessMap;
+ ExcessMap* _excess;
+
+ Tolerance _tolerance;
+
+ bool _phase;
+
+ // constant for treesize
+ static const double _tree_bound = 2;
+ int _max_tree_size;
+
+ //tags for the dynamic tree
+ DynTree *_dt;
+ //datastructure of dyanamic tree
+ IntNodeMap *_dt_index;
+ //datastrucure for solution of the communication between the two class
+ EdgeNodeMap *_dt_edges;
+ //nodeMap for storing the outgoing edge from the node in the tree
+
+ //max of the type Value
+ const Value _max_value;
+
+ public:
+
+ typedef GoldbergTarjan Create;
+
+ ///\name Named template parameters
+
+ ///@{
+
+ template <typename _FlowMap>
+ struct DefFlowMapTraits : public Traits {
+ typedef _FlowMap FlowMap;
+ static FlowMap *createFlowMap(const Graph&) {
+ throw UninitializedParameter();
+ }
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting
+ /// FlowMap type
+ ///
+ /// \ref named-templ-param "Named parameter" for setting FlowMap
+ /// type
+ template <typename _FlowMap>
+ struct DefFlowMap
+ : public GoldbergTarjan<Graph, CapacityMap,
+ DefFlowMapTraits<_FlowMap> > {
+ typedef GoldbergTarjan<Graph, CapacityMap,
+ DefFlowMapTraits<_FlowMap> > Create;
+ };
+
+ template <typename _Elevator>
+ struct DefElevatorTraits : public Traits {
+ typedef _Elevator Elevator;
+ static Elevator *createElevator(const Graph&, int) {
+ throw UninitializedParameter();
+ }
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting
+ /// Elevator type
+ ///
+ /// \ref named-templ-param "Named parameter" for setting Elevator
+ /// type
+ template <typename _Elevator>
+ struct DefElevator
+ : public GoldbergTarjan<Graph, CapacityMap,
+ DefElevatorTraits<_Elevator> > {
+ typedef GoldbergTarjan<Graph, CapacityMap,
+ DefElevatorTraits<_Elevator> > Create;
+ };
+
+ template <typename _Elevator>
+ struct DefStandardElevatorTraits : public Traits {
+ typedef _Elevator Elevator;
+ static Elevator *createElevator(const Graph& graph, int max_level) {
+ return new Elevator(graph, max_level);
+ }
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting
+ /// Elevator type
+ ///
+ /// \ref named-templ-param "Named parameter" for setting Elevator
+ /// type. The Elevator should be standard constructor interface, ie.
+ /// the graph and the maximum level should be passed to it.
+ template <typename _Elevator>
+ struct DefStandardElevator
+ : public GoldbergTarjan<Graph, CapacityMap,
+ DefStandardElevatorTraits<_Elevator> > {
+ typedef GoldbergTarjan<Graph, CapacityMap,
+ DefStandardElevatorTraits<_Elevator> > Create;
+ };
+
+
+ ///\ref Exception for the case when s=t.
+
+ ///\ref Exception for the case when the source equals the target.
+ class InvalidArgument : public lemon::LogicError {
+ public:
+ virtual const char* what() const throw() {
+ return "lemon::GoldbergTarjan::InvalidArgument";
+ }
+ };
+
+ public:
+
+ /// \brief The constructor of the class.
+ ///
+ /// The constructor of the class.
+ /// \param graph The directed graph the algorithm runs on.
+ /// \param capacity The capacity of the edges.
+ /// \param source The source node.
+ /// \param target The target node.
+ /// Except the graph, all of these parameters can be reset by
+ /// calling \ref source, \ref target, \ref capacityMap and \ref
+ /// flowMap, resp.
+ GoldbergTarjan(const Graph& graph, const CapacityMap& capacity,
+ Node source, Node target)
+ : _graph(graph), _capacity(&capacity),
+ _node_num(), _source(source), _target(target),
+ _flow(0), _local_flow(false),
+ _level(0), _local_level(false),
+ _excess(0), _tolerance(),
+ _phase(), _max_tree_size(),
+ _dt(0), _dt_index(0), _dt_edges(0),
+ _max_value(std::numeric_limits<Value>::max()) {
+ if (_source == _target) throw InvalidArgument();
+ }
+
+ /// \brief Destrcutor.
+ ///
+ /// Destructor.
+ ~GoldbergTarjan() {
+ destroyStructures();
+ }
+
+ /// \brief Sets the capacity map.
+ ///
+ /// Sets the capacity map.
+ /// \return \c (*this)
+ GoldbergTarjan& capacityMap(const CapacityMap& map) {
+ _capacity = ↦
+ return *this;
+ }
+
+ /// \brief Sets the flow map.
+ ///
+ /// Sets the flow map.
+ /// \return \c (*this)
+ GoldbergTarjan& flowMap(FlowMap& map) {
+ if (_local_flow) {
+ delete _flow;
+ _local_flow = false;
+ }
+ _flow = ↦
+ return *this;
+ }
+
+ /// \brief Returns the flow map.
+ ///
+ /// \return The flow map.
+ const FlowMap& flowMap() {
+ return *_flow;
+ }
+
+ /// \brief Sets the elevator.
+ ///
+ /// Sets the elevator.
+ /// \return \c (*this)
+ GoldbergTarjan& elevator(Elevator& elevator) {
+ if (_local_level) {
+ delete _level;
+ _local_level = false;
+ }
+ _level = &elevator;
+ return *this;
+ }
+
+ /// \brief Returns the elevator.
+ ///
+ /// \return The elevator.
+ const Elevator& elevator() {
+ return *_level;
+ }
+
+ /// \brief Sets the source node.
+ ///
+ /// Sets the source node.
+ /// \return \c (*this)
+ GoldbergTarjan& source(const Node& node) {
+ _source = node;
+ return *this;
+ }
+
+ /// \brief Sets the target node.
+ ///
+ /// Sets the target node.
+ /// \return \c (*this)
+ GoldbergTarjan& target(const Node& node) {
+ _target = node;
+ return *this;
+ }
+
+ /// \brief Sets the tolerance used by algorithm.
+ ///
+ /// Sets the tolerance used by algorithm.
+ GoldbergTarjan& tolerance(const Tolerance& tolerance) const {
+ _tolerance = tolerance;
+ if (_dt) {
+ _dt->tolerance(_tolerance);
+ }
+ return *this;
+ }
+
+ /// \brief Returns the tolerance used by algorithm.
+ ///
+ /// Returns the tolerance used by algorithm.
+ const Tolerance& tolerance() const {
+ return tolerance;
+ }
+
+
+ private:
+
+ void createStructures() {
+ _node_num = countNodes(_graph);
+
+ _max_tree_size = (double(_node_num) * double(_node_num)) /
+ double(countEdges(_graph));
+
+ if (!_flow) {
+ _flow = Traits::createFlowMap(_graph);
+ _local_flow = true;
+ }
+ if (!_level) {
+ _level = Traits::createElevator(_graph, _node_num);
+ _local_level = true;
+ }
+ if (!_excess) {
+ _excess = new ExcessMap(_graph);
+ }
+ if (!_dt_index && !_dt) {
+ _dt_index = new IntNodeMap(_graph);
+ _dt = new DynTree(*_dt_index, _tolerance);
+ }
+ if (!_dt_edges) {
+ _dt_edges = new EdgeNodeMap(_graph);
+ }
+ }
+
+ void destroyStructures() {
+ if (_local_flow) {
+ delete _flow;
+ }
+ if (_local_level) {
+ delete _level;
+ }
+ if (_excess) {
+ delete _excess;
+ }
+ if (_dt) {
+ delete _dt;
+ }
+ if (_dt_index) {
+ delete _dt_index;
+ }
+ if (_dt_edges) {
+ delete _dt_edges;
+ }
+ }
+
+ bool sendOut(Node n, Value& excess) {
+
+ Node w = _dt->findRoot(n);
+
+ while (w != n) {
+
+ Value rem;
+ Node u = _dt->findCost(n, rem);
+
+ if (_tolerance.positive(rem) && !_level->active(w) && w != _target) {
+ _level->activate(w);
+ }
+
+ if (!_tolerance.less(rem, excess)) {
+ _dt->addCost(n, - excess);
+ _excess->set(w, (*_excess)[w] + excess);
+ excess = 0;
+ return true;
+ }
+
+
+ _dt->addCost(n, - rem);
+
+ excess -= rem;
+ _excess->set(w, (*_excess)[w] + rem);
+
+ _dt->cut(u);
+ _dt->addCost(u, _max_value);
+
+ Edge e = (*_dt_edges)[u];
+ _dt_edges->set(u, INVALID);
+
+ if (u == _graph.source(e)) {
+ _flow->set(e, (*_capacity)[e]);
+ } else {
+ _flow->set(e, 0);
+ }
+
+ w = _dt->findRoot(n);
+ }
+ return false;
+ }
+
+ bool sendIn(Node n, Value& excess) {
+
+ Node w = _dt->findRoot(n);
+
+ while (w != n) {
+
+ Value rem;
+ Node u = _dt->findCost(n, rem);
+
+ if (_tolerance.positive(rem) && !_level->active(w) && w != _source) {
+ _level->activate(w);
+ }
+
+ if (!_tolerance.less(rem, excess)) {
+ _dt->addCost(n, - excess);
+ _excess->set(w, (*_excess)[w] + excess);
+ excess = 0;
+ return true;
+ }
+
+
+ _dt->addCost(n, - rem);
+
+ excess -= rem;
+ _excess->set(w, (*_excess)[w] + rem);
+
+ _dt->cut(u);
+ _dt->addCost(u, _max_value);
+
+ Edge e = (*_dt_edges)[u];
+ _dt_edges->set(u, INVALID);
+
+ if (u == _graph.source(e)) {
+ _flow->set(e, (*_capacity)[e]);
+ } else {
+ _flow->set(e, 0);
+ }
+
+ w = _dt->findRoot(n);
+ }
+ return false;
+ }
+
+ void cutChildren(Node n) {
+
+ for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+
+ Node v = _graph.target(e);
+
+ if ((*_dt_edges)[v] != INVALID && (*_dt_edges)[v] == e) {
+ _dt->cut(v);
+ _dt_edges->set(v, INVALID);
+ Value rem;
+ _dt->findCost(v, rem);
+ _dt->addCost(v, - rem);
+ _dt->addCost(v, _max_value);
+ _flow->set(e, rem);
+
+ }
+ }
+
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+
+ Node v = _graph.source(e);
+
+ if ((*_dt_edges)[v] != INVALID && (*_dt_edges)[v] == e) {
+ _dt->cut(v);
+ _dt_edges->set(v, INVALID);
+ Value rem;
+ _dt->findCost(v, rem);
+ _dt->addCost(v, - rem);
+ _dt->addCost(v, _max_value);
+ _flow->set(e, (*_capacity)[e] - rem);
+
+ }
+ }
+ }
+
+ void extractTrees() {
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+
+ Node w = _dt->findRoot(n);
+
+ while (w != n) {
+
+ Value rem;
+ Node u = _dt->findCost(n, rem);
+
+ _dt->cut(u);
+ _dt->addCost(u, - rem);
+ _dt->addCost(u, _max_value);
+
+ Edge e = (*_dt_edges)[u];
+ _dt_edges->set(u, INVALID);
+
+ if (u == _graph.source(e)) {
+ _flow->set(e, (*_capacity)[e] - rem);
+ } else {
+ _flow->set(e, rem);
+ }
+
+ w = _dt->findRoot(n);
+ }
+ }
+ }
+
+ public:
+
+ /// \name Execution control The simplest way to execute the
+ /// algorithm is to use one of the member functions called \c
+ /// run().
+ /// \n
+ /// If you need more control on initial solution or
+ /// execution then you have to call one \ref init() function and then
+ /// the startFirstPhase() and if you need the startSecondPhase().
+
+ ///@{
+
+ /// \brief Initializes the internal data structures.
+ ///
+ /// Initializes the internal data structures.
+ ///
+ void init() {
+ createStructures();
+
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ _excess->set(n, 0);
+ }
+
+ for (EdgeIt e(_graph); e != INVALID; ++e) {
+ _flow->set(e, 0);
+ }
+
+ _dt->clear();
+ for (NodeIt v(_graph); v != INVALID; ++v) {
+ (*_dt_edges)[v] = INVALID;
+ _dt->makeTree(v);
+ _dt->addCost(v, _max_value);
+ }
+
+ typename Graph::template NodeMap<bool> reached(_graph, false);
+
+ _level->initStart();
+ _level->initAddItem(_target);
+
+ std::vector<Node> queue;
+ reached.set(_source, true);
+
+ queue.push_back(_target);
+ reached.set(_target, true);
+ while (!queue.empty()) {
+ _level->initNewLevel();
+ std::vector<Node> nqueue;
+ for (int i = 0; i < int(queue.size()); ++i) {
+ Node n = queue[i];
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Node u = _graph.source(e);
+ if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
+ reached.set(u, true);
+ _level->initAddItem(u);
+ nqueue.push_back(u);
+ }
+ }
+ }
+ queue.swap(nqueue);
+ }
+ _level->initFinish();
+
+ for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
+ if (_tolerance.positive((*_capacity)[e])) {
+ Node u = _graph.target(e);
+ if ((*_level)[u] == _level->maxLevel()) continue;
+ _flow->set(e, (*_capacity)[e]);
+ _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
+ if (u != _target && !_level->active(u)) {
+ _level->activate(u);
+ }
+ }
+ }
+ }
+
+ /// \brief Starts the first phase of the preflow algorithm.
+ ///
+ /// The preflow algorithm consists of two phases, this method runs
+ /// the first phase. After the first phase the maximum flow value
+ /// and a minimum value cut can already be computed, although a
+ /// maximum flow is not yet obtained. So after calling this method
+ /// \ref flowValue() returns the value of a maximum flow and \ref
+ /// minCut() returns a minimum cut.
+ /// \pre One of the \ref init() functions should be called.
+ void startFirstPhase() {
+ _phase = true;
+ Node n;
+
+ while ((n = _level->highestActive()) != INVALID) {
+ Value excess = (*_excess)[n];
+ int level = _level->highestActiveLevel();
+ int new_level = _level->maxLevel();
+
+ if (_dt->findRoot(n) != n) {
+ if (sendOut(n, excess)) goto no_more_push;
+ }
+
+ for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_capacity)[e] - (*_flow)[e];
+ Node v = _graph.target(e);
+
+ if (!_tolerance.positive(rem) && (*_dt_edges)[v] != e) continue;
+
+ if ((*_level)[v] < level) {
+
+ if (_dt->findSize(n) + _dt->findSize(v) <
+ _tree_bound * _max_tree_size) {
+ _dt->addCost(n, -_max_value);
+ _dt->addCost(n, rem);
+ _dt->link(n, v);
+ _dt_edges->set(n, e);
+ if (sendOut(n, excess)) goto no_more_push;
+ } else {
+ if (!_level->active(v) && v != _target) {
+ _level->activate(v);
+ }
+ if (!_tolerance.less(rem, excess)) {
+ _flow->set(e, (*_flow)[e] + excess);
+ _excess->set(v, (*_excess)[v] + excess);
+ excess = 0;
+ goto no_more_push;
+ } else {
+ excess -= rem;
+ _excess->set(v, (*_excess)[v] + rem);
+ _flow->set(e, (*_capacity)[e]);
+ }
+ }
+ } else if (new_level > (*_level)[v]) {
+ new_level = (*_level)[v];
+ }
+ }
+
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_flow)[e];
+ Node v = _graph.source(e);
+ if (!_tolerance.positive(rem) && (*_dt_edges)[v] != e) continue;
+
+ if ((*_level)[v] < level) {
+
+ if (_dt->findSize(n) + _dt->findSize(v) <
+ _tree_bound * _max_tree_size) {
+ _dt->addCost(n, - _max_value);
+ _dt->addCost(n, rem);
+ _dt->link(n, v);
+ _dt_edges->set(n, e);
+ if (sendOut(n, excess)) goto no_more_push;
+ } else {
+ if (!_level->active(v) && v != _target) {
+ _level->activate(v);
+ }
+ if (!_tolerance.less(rem, excess)) {
+ _flow->set(e, (*_flow)[e] - excess);
+ _excess->set(v, (*_excess)[v] + excess);
+ excess = 0;
+ goto no_more_push;
+ } else {
+ excess -= rem;
+ _excess->set(v, (*_excess)[v] + rem);
+ _flow->set(e, 0);
+ }
+ }
+ } else if (new_level > (*_level)[v]) {
+ new_level = (*_level)[v];
+ }
+ }
+
+ no_more_push:
+
+ _excess->set(n, excess);
+
+ if (excess != 0) {
+ cutChildren(n);
+ if (new_level + 1 < _level->maxLevel()) {
+ _level->liftHighestActive(new_level + 1);
+ } else {
+ _level->liftHighestActiveToTop();
+ }
+ if (_level->emptyLevel(level)) {
+ _level->liftToTop(level);
+ }
+ } else {
+ _level->deactivate(n);
+ }
+ }
+ extractTrees();
+ }
+
+ /// \brief Starts the second phase of the preflow algorithm.
+ ///
+ /// The preflow algorithm consists of two phases, this method runs
+ /// the second phase. After calling \ref init() and \ref
+ /// startFirstPhase() and then \ref startSecondPhase(), \ref
+ /// flowMap() return a maximum flow, \ref flowValue() returns the
+ /// value of a maximum flow, \ref minCut() returns a minimum cut
+ /// \pre The \ref init() and startFirstPhase() functions should be
+ /// called before.
+ void startSecondPhase() {
+ _phase = false;
+
+ typename Graph::template NodeMap<bool> reached(_graph, false);
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ reached.set(n, (*_level)[n] < _level->maxLevel());
+ }
+
+ _level->initStart();
+ _level->initAddItem(_source);
+
+ std::vector<Node> queue;
+ queue.push_back(_source);
+ reached.set(_source, true);
+
+ while (!queue.empty()) {
+ _level->initNewLevel();
+ std::vector<Node> nqueue;
+ for (int i = 0; i < int(queue.size()); ++i) {
+ Node n = queue[i];
+ for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Node v = _graph.target(e);
+ if (!reached[v] && _tolerance.positive((*_flow)[e])) {
+ reached.set(v, true);
+ _level->initAddItem(v);
+ nqueue.push_back(v);
+ }
+ }
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Node u = _graph.source(e);
+ if (!reached[u] &&
+ _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
+ reached.set(u, true);
+ _level->initAddItem(u);
+ nqueue.push_back(u);
+ }
+ }
+ }
+ queue.swap(nqueue);
+ }
+ _level->initFinish();
+
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ if ((*_excess)[n] > 0 && _target != n) {
+ _level->activate(n);
+ }
+ }
+
+ Node n;
+
+ while ((n = _level->highestActive()) != INVALID) {
+ Value excess = (*_excess)[n];
+ int level = _level->highestActiveLevel();
+ int new_level = _level->maxLevel();
+
+ if (_dt->findRoot(n) != n) {
+ if (sendIn(n, excess)) goto no_more_push;
+ }
+
+ for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_capacity)[e] - (*_flow)[e];
+ Node v = _graph.target(e);
+
+ if (!_tolerance.positive(rem) && (*_dt_edges)[v] != e) continue;
+
+ if ((*_level)[v] < level) {
+
+ if (_dt->findSize(n) + _dt->findSize(v) <
+ _tree_bound * _max_tree_size) {
+ _dt->addCost(n, -_max_value);
+ _dt->addCost(n, rem);
+ _dt->link(n, v);
+ _dt_edges->set(n, e);
+ if (sendIn(n, excess)) goto no_more_push;
+ } else {
+ if (!_level->active(v) && v != _source) {
+ _level->activate(v);
+ }
+ if (!_tolerance.less(rem, excess)) {
+ _flow->set(e, (*_flow)[e] + excess);
+ _excess->set(v, (*_excess)[v] + excess);
+ excess = 0;
+ goto no_more_push;
+ } else {
+ excess -= rem;
+ _excess->set(v, (*_excess)[v] + rem);
+ _flow->set(e, (*_capacity)[e]);
+ }
+ }
+ } else if (new_level > (*_level)[v]) {
+ new_level = (*_level)[v];
+ }
+ }
+
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_flow)[e];
+ Node v = _graph.source(e);
+ if (!_tolerance.positive(rem) && (*_dt_edges)[v] != e) continue;
+
+ if ((*_level)[v] < level) {
+
+ if (_dt->findSize(n) + _dt->findSize(v) <
+ _tree_bound * _max_tree_size) {
+ _dt->addCost(n, - _max_value);
+ _dt->addCost(n, rem);
+ _dt->link(n, v);
+ _dt_edges->set(n, e);
+ if (sendIn(n, excess)) goto no_more_push;
+ } else {
+ if (!_level->active(v) && v != _source) {
+ _level->activate(v);
+ }
+ if (!_tolerance.less(rem, excess)) {
+ _flow->set(e, (*_flow)[e] - excess);
+ _excess->set(v, (*_excess)[v] + excess);
+ excess = 0;
+ goto no_more_push;
+ } else {
+ excess -= rem;
+ _excess->set(v, (*_excess)[v] + rem);
+ _flow->set(e, 0);
+ }
+ }
+ } else if (new_level > (*_level)[v]) {
+ new_level = (*_level)[v];
+ }
+ }
+
+ no_more_push:
+
+ _excess->set(n, excess);
+
+ if (excess != 0) {
+ cutChildren(n);
+ if (new_level + 1 < _level->maxLevel()) {
+ _level->liftHighestActive(new_level + 1);
+ } else {
+ _level->liftHighestActiveToTop();
+ }
+ if (_level->emptyLevel(level)) {
+ _level->liftToTop(level);
+ }
+ } else {
+ _level->deactivate(n);
+ }
+ }
+ extractTrees();
+ }
+
+ /// \brief Runs the Goldberg-Tarjan algorithm.
+ ///
+ /// Runs the Goldberg-Tarjan algorithm.
+ /// \note pf.run() is just a shortcut of the following code.
+ /// \code
+ /// pf.init();
+ /// pf.startFirstPhase();
+ /// pf.startSecondPhase();
+ /// \endcode
+ void run() {
+ init();
+ startFirstPhase();
+ startSecondPhase();
+ }
+
+ /// \brief Runs the Goldberg-Tarjan algorithm to compute the minimum cut.
+ ///
+ /// Runs the Goldberg-Tarjan algorithm to compute the minimum cut.
+ /// \note pf.runMinCut() is just a shortcut of the following code.
+ /// \code
+ /// pf.init();
+ /// pf.startFirstPhase();
+ /// \endcode
+ void runMinCut() {
+ init();
+ startFirstPhase();
+ }
+
+ /// @}
+
+ /// \name Query Functions
+ /// The result of the %Dijkstra algorithm can be obtained using these
+ /// functions.\n
+ /// Before the use of these functions,
+ /// either run() or start() must be called.
+
+ ///@{
+
+ /// \brief Returns the value of the maximum flow.
+ ///
+ /// Returns the value of the maximum flow by returning the excess
+ /// of the target node \c t. This value equals to the value of
+ /// the maximum flow already after the first phase.
+ Value flowValue() const {
+ return (*_excess)[_target];
+ }
+
+ /// \brief Returns true when the node is on the source side of minimum cut.
+ ///
+ /// Returns true when the node is on the source side of minimum
+ /// cut. This method can be called both after running \ref
+ /// startFirstPhase() and \ref startSecondPhase().
+ bool minCut(const Node& node) const {
+ return ((*_level)[node] == _level->maxLevel()) == _phase;
+ }
+
+ /// \brief Returns a minimum value cut.
+ ///
+ /// Sets the \c cutMap to the characteristic vector of a minimum value
+ /// cut. This method can be called both after running \ref
+ /// startFirstPhase() and \ref startSecondPhase(). The result after second
+ /// phase could be changed slightly if inexact computation is used.
+ /// \pre The \c cutMap should be a bool-valued node-map.
+ template <typename CutMap>
+ void minCutMap(CutMap& cutMap) const {
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ cutMap.set(n, minCut(n));
+ }
+ }
+
+ /// \brief Returns the flow on the edge.
+ ///
+ /// Sets the \c flowMap to the flow on the edges. This method can
+ /// be called after the second phase of algorithm.
+ Value flow(const Edge& edge) const {
+ return (*_flow)[edge];
+ }
+
+ /// @}
+
+ };
+
+} //namespace lemon
+
+#endif
Modified: lemon/trunk/lemon/preflow.h
==============================================================================
--- lemon/trunk/lemon/preflow.h (original)
+++ lemon/trunk/lemon/preflow.h Sat Nov 17 21:58:11 2007
@@ -19,897 +19,897 @@
#ifndef LEMON_PREFLOW_H
#define LEMON_PREFLOW_H
-#include <vector>
-#include <queue>
-
#include <lemon/error.h>
#include <lemon/bits/invalid.h>
#include <lemon/tolerance.h>
#include <lemon/maps.h>
#include <lemon/graph_utils.h>
+#include <lemon/elevator.h>
/// \file
/// \ingroup max_flow
/// \brief Implementation of the preflow algorithm.
-namespace lemon {
-
- ///\ingroup max_flow
- ///\brief %Preflow algorithms class.
- ///
- ///This class provides an implementation of the \e preflow \e
- ///algorithm producing a flow of maximum value in a directed
- ///graph. The preflow algorithms are the fastest known max flow algorithms.
- ///The \e source node, the \e target node, the \e
- ///capacity of the edges and the \e starting \e flow value of the
- ///edges should be passed to the algorithm through the
- ///constructor. It is possible to change these quantities using the
- ///functions \ref source, \ref target, \ref capacityMap and \ref
- ///flowMap.
- ///
- ///After running \ref lemon::Preflow::phase1() "phase1()"
- ///or \ref lemon::Preflow::run() "run()", the maximal flow
- ///value can be obtained by calling \ref flowValue(). The minimum
- ///value cut can be written into a <tt>bool</tt> node map by
- ///calling \ref minCut(). (\ref minMinCut() and \ref maxMinCut() writes
- ///the inclusionwise minimum and maximum of the minimum value cuts,
- ///resp.)
- ///
- ///\param Graph The directed graph type the algorithm runs on.
- ///\param Num The number type of the capacities and the flow values.
- ///\param CapacityMap The capacity map type.
- ///\param FlowMap The flow map type.
- ///\param Tol The tolerance type.
+namespace lemon {
+
+ /// \brief Default traits class of Preflow class.
///
- ///\author Jacint Szabo
- ///\todo Second template parameter is superfluous
- template <typename Graph, typename Num,
- typename CapacityMap=typename Graph::template EdgeMap<Num>,
- typename FlowMap=typename Graph::template EdgeMap<Num>,
- typename Tol=Tolerance<Num> >
- class Preflow {
- protected:
- typedef typename Graph::Node Node;
- typedef typename Graph::NodeIt NodeIt;
- typedef typename Graph::EdgeIt EdgeIt;
- typedef typename Graph::OutEdgeIt OutEdgeIt;
- typedef typename Graph::InEdgeIt InEdgeIt;
-
- typedef typename Graph::template NodeMap<Node> NNMap;
- typedef typename std::vector<Node> VecNode;
-
- const Graph* _g;
- Node _source;
- Node _target;
- const CapacityMap* _capacity;
- FlowMap* _flow;
+ /// Default traits class of Preflow class.
+ /// \param _Graph Graph type.
+ /// \param _CapacityMap Type of capacity map.
+ template <typename _Graph, typename _CapacityMap>
+ struct PreflowDefaultTraits {
- Tol _surely;
-
- int _node_num; //the number of nodes of G
-
- typename Graph::template NodeMap<int> level;
- typename Graph::template NodeMap<Num> excess;
-
- // constants used for heuristics
- static const int H0=20;
- static const int H1=1;
-
- public:
+ /// \brief The graph type the algorithm runs on.
+ typedef _Graph Graph;
- ///\ref Exception for the case when s=t.
-
- ///\ref Exception for the case when the source equals the target.
- class InvalidArgument : public lemon::LogicError {
- public:
- virtual const char* what() const throw() {
- return "lemon::Preflow::InvalidArgument";
- }
- };
-
-
- ///Indicates the property of the starting flow map.
-
- ///Indicates the property of the starting flow map.
+ /// \brief The type of the map that stores the edge capacities.
///
- enum FlowEnum{
- ///indicates an unspecified edge map. \c flow will be
- ///set to the constant zero flow in the beginning of
- ///the algorithm in this case.
- NO_FLOW,
- ///constant zero flow
- ZERO_FLOW,
- ///any flow, i.e. the sum of the in-flows equals to
- ///the sum of the out-flows in every node except the \c source and
- ///the \c target.
- GEN_FLOW,
- ///any preflow, i.e. the sum of the in-flows is at
- ///least the sum of the out-flows in every node except the \c source.
- PRE_FLOW
- };
+ /// The type of the map that stores the edge capacities.
+ /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
+ typedef _CapacityMap CapacityMap;
- ///Indicates the state of the preflow algorithm.
+ /// \brief The type of the length of the edges.
+ typedef typename CapacityMap::Value Value;
- ///Indicates the state of the preflow algorithm.
+ /// \brief The map type that stores the flow values.
///
- enum StatusEnum {
- ///before running the algorithm or
- ///at an unspecified state.
- AFTER_NOTHING,
- ///right after running \ref phase1()
- AFTER_PREFLOW_PHASE_1,
- ///after running \ref phase2()
- AFTER_PREFLOW_PHASE_2
- };
-
- protected:
- FlowEnum flow_prop;
- StatusEnum status; // Do not needle this flag only if necessary.
-
- public:
- ///The constructor of the class.
+ /// The map type that stores the flow values.
+ /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
+ typedef typename Graph::template EdgeMap<Value> FlowMap;
- ///The constructor of the class.
- ///\param _gr The directed graph the algorithm runs on.
- ///\param _s The source node.
- ///\param _t The target node.
- ///\param _cap The capacity of the edges.
- ///\param _f The flow of the edges.
- ///\param _sr Tol class.
- ///Except the graph, all of these parameters can be reset by
- ///calling \ref source, \ref target, \ref capacityMap and \ref
- ///flowMap, resp.
- Preflow(const Graph& _gr, Node _s, Node _t,
- const CapacityMap& _cap, FlowMap& _f,
- const Tol &_sr=Tol()) :
- _g(&_gr), _source(_s), _target(_t), _capacity(&_cap),
- _flow(&_f), _surely(_sr),
- _node_num(countNodes(_gr)), level(_gr), excess(_gr,0),
- flow_prop(NO_FLOW), status(AFTER_NOTHING) {
- if ( _source==_target )
- throw InvalidArgument();
+ /// \brief Instantiates a FlowMap.
+ ///
+ /// This function instantiates a \ref FlowMap.
+ /// \param graph The graph, to which we would like to define the flow map.
+ static FlowMap* createFlowMap(const Graph& graph) {
+ return new FlowMap(graph);
}
-
- ///Give a reference to the tolerance handler class
-
- ///Give a reference to the tolerance handler class
- ///\sa Tolerance
- Tol &tolerance() { return _surely; }
- ///Runs the preflow algorithm.
-
- ///Runs the preflow algorithm.
+ /// \brief The eleavator type used by Preflow algorithm.
+ ///
+ /// The elevator type used by Preflow algorithm.
///
- void run() {
- phase1(flow_prop);
- phase2();
- }
+ /// \sa Elevator
+ /// \sa LinkedElevator
+ typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
- ///Runs the preflow algorithm.
-
- ///Runs the preflow algorithm.
- ///\pre The starting flow map must be
- /// - a constant zero flow if \c fp is \c ZERO_FLOW,
- /// - an arbitrary flow if \c fp is \c GEN_FLOW,
- /// - an arbitrary preflow if \c fp is \c PRE_FLOW,
- /// - any map if \c fp is NO_FLOW.
- ///If the starting flow map is a flow or a preflow then
- ///the algorithm terminates faster.
- void run(FlowEnum fp) {
- flow_prop=fp;
- run();
+ /// \brief Instantiates an Elevator.
+ ///
+ /// This function instantiates a \ref Elevator.
+ /// \param graph The graph, to which we would like to define the elevator.
+ /// \param max_level The maximum level of the elevator.
+ static Elevator* createElevator(const Graph& graph, int max_level) {
+ return new Elevator(graph, max_level);
}
-
- ///Runs the first phase of the preflow algorithm.
- ///The preflow algorithm consists of two phases, this method runs
- ///the first phase. After the first phase the maximum flow value
- ///and a minimum value cut can already be computed, although a
- ///maximum flow is not yet obtained. So after calling this method
- ///\ref flowValue returns the value of a maximum flow and \ref
- ///minCut returns a minimum cut.
- ///\warning \ref minMinCut and \ref maxMinCut do not give minimum
- ///value cuts unless calling \ref phase2.
- ///\warning A real flow map (i.e. not \ref lemon::NullMap "NullMap")
- ///is needed for this phase.
- ///\pre The starting flow must be
- ///- a constant zero flow if \c fp is \c ZERO_FLOW,
- ///- an arbitary flow if \c fp is \c GEN_FLOW,
- ///- an arbitary preflow if \c fp is \c PRE_FLOW,
- ///- any map if \c fp is NO_FLOW.
- void phase1(FlowEnum fp)
- {
- flow_prop=fp;
- phase1();
- }
+ /// \brief The tolerance used by the algorithm
+ ///
+ /// The tolerance used by the algorithm to handle inexact computation.
+ typedef Tolerance<Value> Tolerance;
+ };
+
+
+ /// \ingroup max_flow
+ ///
+ /// \brief %Preflow algorithms class.
+ ///
+ /// This class provides an implementation of the Goldberg's \e
+ /// preflow \e algorithm producing a flow of maximum value in a
+ /// directed graph. The preflow algorithms are the fastest known max
+ /// flow algorithms. The current implementation use a mixture of the
+ /// \e "highest label" and the \e "bound decrease" heuristics.
+ /// The worst case time complexity of the algorithm is \f$O(n^3)\f$.
+ ///
+ /// The algorithm consists from two phases. After the first phase
+ /// the maximal flow value and the minimum cut can be obtained. The
+ /// second phase constructs the feasible maximum flow on each edge.
+ ///
+ /// \param _Graph The directed graph type the algorithm runs on.
+ /// \param _CapacityMap The flow map type.
+ /// \param _Traits Traits class to set various data types used by
+ /// the algorithm. The default traits class is \ref
+ /// PreflowDefaultTraits. See \ref PreflowDefaultTraits for the
+ /// documentation of a %Preflow traits class.
+ ///
+ ///\author Jacint Szabo and Balazs Dezso
+#ifdef DOXYGEN
+ template <typename _Graph, typename _CapacityMap, typename _Traits>
+#else
+ template <typename _Graph,
+ typename _CapacityMap = typename _Graph::template EdgeMap<int>,
+ typename _Traits = PreflowDefaultTraits<_Graph, _CapacityMap> >
+#endif
+ class Preflow {
+ public:
- ///Runs the first phase of the preflow algorithm.
+ typedef _Traits Traits;
+ typedef typename Traits::Graph Graph;
+ typedef typename Traits::CapacityMap CapacityMap;
+ typedef typename Traits::Value Value;
+
+ typedef typename Traits::FlowMap FlowMap;
+ typedef typename Traits::Elevator Elevator;
+ typedef typename Traits::Tolerance Tolerance;
- ///The preflow algorithm consists of two phases, this method runs
- ///the first phase. After the first phase the maximum flow value
- ///and a minimum value cut can already be computed, although a
- ///maximum flow is not yet obtained. So after calling this method
- ///\ref flowValue returns the value of a maximum flow and \ref
- ///minCut returns a minimum cut.
- ///\warning \ref minMinCut() and \ref maxMinCut() do not
- ///give minimum value cuts unless calling \ref phase2().
- ///\warning A real flow map (i.e. not \ref lemon::NullMap "NullMap")
- ///is needed for this phase.
- void phase1()
- {
- int heur0=int(H0*_node_num); //time while running 'bound decrease'
- int heur1=int(H1*_node_num); //time while running 'highest label'
- int heur=heur1; //starting time interval (#of relabels)
- int numrelabel=0;
-
- bool what_heur=1;
- //It is 0 in case 'bound decrease' and 1 in case 'highest label'
-
- bool end=false;
- //Needed for 'bound decrease', true means no active
- //nodes are above bound b.
-
- int k=_node_num-2; //bound on the highest level under n containing a node
- int b=k; //bound on the highest level under n containing an active node
-
- VecNode first(_node_num, INVALID);
- NNMap next(*_g, INVALID);
-
- NNMap left(*_g, INVALID);
- NNMap right(*_g, INVALID);
- VecNode level_list(_node_num,INVALID);
- //List of the nodes in level i<n, set to n.
-
- preflowPreproc(first, next, level_list, left, right);
-
- //Push/relabel on the highest level active nodes.
- while ( true ) {
- if ( b == 0 ) {
- if ( !what_heur && !end && k > 0 ) {
- b=k;
- end=true;
- } else break;
- }
-
- if ( first[b]==INVALID ) --b;
- else {
- end=false;
- Node w=first[b];
- first[b]=next[w];
- int newlevel=push(w, next, first);
- if ( excess[w] != 0 ) {
- relabel(w, newlevel, first, next, level_list,
- left, right, b, k, what_heur);
- }
-
- ++numrelabel;
- if ( numrelabel >= heur ) {
- numrelabel=0;
- if ( what_heur ) {
- what_heur=0;
- heur=heur0;
- end=false;
- } else {
- what_heur=1;
- heur=heur1;
- b=k;
- }
- }
- }
+ /// \brief \ref Exception for uninitialized parameters.
+ ///
+ /// This error represents problems in the initialization
+ /// of the parameters of the algorithms.
+ class UninitializedParameter : public lemon::UninitializedParameter {
+ public:
+ virtual const char* what() const throw() {
+ return "lemon::Preflow::UninitializedParameter";
}
- flow_prop=PRE_FLOW;
- status=AFTER_PREFLOW_PHASE_1;
- }
- // Heuristics:
- // 2 phase
- // gap
- // list 'level_list' on the nodes on level i implemented by hand
- // stack 'active' on the active nodes on level i
- // runs heuristic 'highest label' for H1*n relabels
- // runs heuristic 'bound decrease' for H0*n relabels,
- // starts with 'highest label'
- // Parameters H0 and H1 are initialized to 20 and 1.
-
-
- ///Runs the second phase of the preflow algorithm.
-
- ///The preflow algorithm consists of two phases, this method runs
- ///the second phase. After calling \ref phase1() and then
- ///\ref phase2(),
- /// \ref flowMap() return a maximum flow, \ref flowValue
- ///returns the value of a maximum flow, \ref minCut returns a
- ///minimum cut, while the methods \ref minMinCut and \ref
- ///maxMinCut return the inclusionwise minimum and maximum cuts of
- ///minimum value, resp. \pre \ref phase1 must be called before.
- ///
- /// \todo The inexact computation can cause positive excess on a set of
- /// unpushable nodes. We may have to watch the empty level in this case
- /// due to avoid the terrible long running time.
- void phase2()
- {
+ };
- int k=_node_num-2; //bound on the highest level under n containing a node
- int b=k; //bound on the highest level under n of an active node
+ private:
-
- VecNode first(_node_num, INVALID);
- NNMap next(*_g, INVALID);
- level.set(_source,0);
- std::queue<Node> bfs_queue;
- bfs_queue.push(_source);
+ GRAPH_TYPEDEFS(typename Graph);
- while ( !bfs_queue.empty() ) {
+ const Graph& _graph;
+ const CapacityMap* _capacity;
- Node v=bfs_queue.front();
- bfs_queue.pop();
- int l=level[v]+1;
+ int _node_num;
- for(InEdgeIt e(*_g,v); e!=INVALID; ++e) {
- if ( !_surely.positive((*_capacity)[e] - (*_flow)[e])) continue;
- Node u=_g->source(e);
- if ( level[u] >= _node_num ) {
- bfs_queue.push(u);
- level.set(u, l);
- if ( excess[u] != 0 ) {
- next.set(u,first[l]);
- first[l]=u;
- }
- }
- }
+ Node _source, _target;
- for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) {
- if ( !_surely.positive((*_flow)[e]) ) continue;
- Node u=_g->target(e);
- if ( level[u] >= _node_num ) {
- bfs_queue.push(u);
- level.set(u, l);
- if ( excess[u] != 0 ) {
- next.set(u,first[l]);
- first[l]=u;
- }
- }
- }
- }
- b=_node_num-2;
+ FlowMap* _flow;
+ bool _local_flow;
- while ( true ) {
+ Elevator* _level;
+ bool _local_level;
- if ( b == 0 ) break;
- if ( first[b]==INVALID ) --b;
- else {
- Node w=first[b];
- first[b]=next[w];
- int newlevel=push(w,next, first);
-
- if ( newlevel == _node_num) {
- excess.set(w, 0);
- level.set(w,_node_num);
- }
- //relabel
- if ( excess[w] != 0 ) {
- level.set(w,++newlevel);
- next.set(w,first[newlevel]);
- first[newlevel]=w;
- b=newlevel;
- }
- }
- } // while(true)
- flow_prop=GEN_FLOW;
- status=AFTER_PREFLOW_PHASE_2;
- }
+ typedef typename Graph::template NodeMap<Value> ExcessMap;
+ ExcessMap* _excess;
- /// Returns the value of the maximum flow.
+ Tolerance _tolerance;
- /// Returns the value of the maximum flow by returning the excess
- /// of the target node \c t. This value equals to the value of
- /// the maximum flow already after running \ref phase1.
- Num flowValue() const {
- return excess[_target];
- }
+ bool _phase;
+ void createStructures() {
+ _node_num = countNodes(_graph);
- ///Returns a minimum value cut.
-
- ///Sets \c M to the characteristic vector of a minimum value
- ///cut. This method can be called both after running \ref
- ///phase1 and \ref phase2. It is much faster after
- ///\ref phase1. \pre M should be a bool-valued node-map. \pre
- ///If \ref minCut() is called after \ref phase2() then M should
- ///be initialized to false.
- template<typename _CutMap>
- void minCut(_CutMap& M) const {
- switch ( status ) {
- case AFTER_PREFLOW_PHASE_1:
- for(NodeIt v(*_g); v!=INVALID; ++v) {
- if (level[v] < _node_num) {
- M.set(v, false);
- } else {
- M.set(v, true);
- }
- }
- break;
- case AFTER_PREFLOW_PHASE_2:
- minMinCut(M);
- break;
- case AFTER_NOTHING:
- break;
+ if (!_flow) {
+ _flow = Traits::createFlowMap(_graph);
+ _local_flow = true;
+ }
+ if (!_level) {
+ _level = Traits::createElevator(_graph, _node_num);
+ _local_level = true;
+ }
+ if (!_excess) {
+ _excess = new ExcessMap(_graph);
}
}
- ///Returns the inclusionwise minimum of the minimum value cuts.
-
- ///Sets \c M to the characteristic vector of the minimum value cut
- ///which is inclusionwise minimum. It is computed by processing a
- ///bfs from the source node \c s in the residual graph. \pre M
- ///should be a node map of bools initialized to false. \pre \ref
- ///phase2 should already be run.
- template<typename _CutMap>
- void minMinCut(_CutMap& M) const {
-
- std::queue<Node> queue;
- M.set(_source,true);
- queue.push(_source);
-
- while (!queue.empty()) {
- Node w=queue.front();
- queue.pop();
-
- for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
- Node v=_g->target(e);
- if (!M[v] && _surely.positive((*_capacity)[e] -(*_flow)[e]) ) {
- queue.push(v);
- M.set(v, true);
- }
- }
-
- for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
- Node v=_g->source(e);
- if (!M[v] && _surely.positive((*_flow)[e]) ) {
- queue.push(v);
- M.set(v, true);
- }
- }
+ void destroyStructures() {
+ if (_local_flow) {
+ delete _flow;
+ }
+ if (_local_level) {
+ delete _level;
+ }
+ if (_excess) {
+ delete _excess;
}
}
-
- ///Returns the inclusionwise maximum of the minimum value cuts.
- ///Sets \c M to the characteristic vector of the minimum value cut
- ///which is inclusionwise maximum. It is computed by processing a
- ///backward bfs from the target node \c t in the residual graph.
- ///\pre \ref phase2() or run() should already be run.
- template<typename _CutMap>
- void maxMinCut(_CutMap& M) const {
+ public:
- for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true);
+ typedef Preflow Create;
- std::queue<Node> queue;
+ ///\name Named template parameters
- M.set(_target,false);
- queue.push(_target);
+ ///@{
- while (!queue.empty()) {
- Node w=queue.front();
- queue.pop();
+ template <typename _FlowMap>
+ struct DefFlowMapTraits : public Traits {
+ typedef _FlowMap FlowMap;
+ static FlowMap *createFlowMap(const Graph&) {
+ throw UninitializedParameter();
+ }
+ };
- for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
- Node v=_g->source(e);
- if (M[v] && _surely.positive((*_capacity)[e] - (*_flow)[e]) ) {
- queue.push(v);
- M.set(v, false);
- }
- }
+ /// \brief \ref named-templ-param "Named parameter" for setting
+ /// FlowMap type
+ ///
+ /// \ref named-templ-param "Named parameter" for setting FlowMap
+ /// type
+ template <typename _FlowMap>
+ struct DefFlowMap
+ : public Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > {
+ typedef Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > Create;
+ };
- for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
- Node v=_g->target(e);
- if (M[v] && _surely.positive((*_flow)[e]) ) {
- queue.push(v);
- M.set(v, false);
- }
- }
+ template <typename _Elevator>
+ struct DefElevatorTraits : public Traits {
+ typedef _Elevator Elevator;
+ static Elevator *createElevator(const Graph&, int) {
+ throw UninitializedParameter();
}
- }
+ };
- ///Sets the source node to \c _s.
+ /// \brief \ref named-templ-param "Named parameter" for setting
+ /// Elevator type
+ ///
+ /// \ref named-templ-param "Named parameter" for setting Elevator
+ /// type
+ template <typename _Elevator>
+ struct DefElevator
+ : public Preflow<Graph, CapacityMap, DefElevatorTraits<_Elevator> > {
+ typedef Preflow<Graph, CapacityMap, DefElevatorTraits<_Elevator> > Create;
+ };
- ///Sets the source node to \c _s.
- ///
- void source(Node _s) {
- _source=_s;
- if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
- status=AFTER_NOTHING;
- }
+ template <typename _Elevator>
+ struct DefStandardElevatorTraits : public Traits {
+ typedef _Elevator Elevator;
+ static Elevator *createElevator(const Graph& graph, int max_level) {
+ return new Elevator(graph, max_level);
+ }
+ };
- ///Returns the source node.
+ /// \brief \ref named-templ-param "Named parameter" for setting
+ /// Elevator type
+ ///
+ /// \ref named-templ-param "Named parameter" for setting Elevator
+ /// type. The Elevator should be standard constructor interface, ie.
+ /// the graph and the maximum level should be passed to it.
+ template <typename _Elevator>
+ struct DefStandardElevator
+ : public Preflow<Graph, CapacityMap,
+ DefStandardElevatorTraits<_Elevator> > {
+ typedef Preflow<Graph, CapacityMap,
+ DefStandardElevatorTraits<_Elevator> > Create;
+ };
- ///Returns the source node.
- ///
- Node source() const {
- return _source;
- }
+ /// @}
- ///Sets the target node to \c _t.
+ /// \brief The constructor of the class.
+ ///
+ /// The constructor of the class.
+ /// \param graph The directed graph the algorithm runs on.
+ /// \param capacity The capacity of the edges.
+ /// \param source The source node.
+ /// \param target The target node.
+ Preflow(const Graph& graph, const CapacityMap& capacity,
+ Node source, Node target)
+ : _graph(graph), _capacity(&capacity),
+ _node_num(0), _source(source), _target(target),
+ _flow(0), _local_flow(false),
+ _level(0), _local_level(false),
+ _excess(0), _tolerance(), _phase() {}
- ///Sets the target node to \c _t.
+ /// \brief Destrcutor.
///
- void target(Node _t) {
- _target=_t;
- if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
- status=AFTER_NOTHING;
+ /// Destructor.
+ ~Preflow() {
+ destroyStructures();
}
- ///Returns the target node.
-
- ///Returns the target node.
- ///
- Node target() const {
- return _target;
+ /// \brief Sets the capacity map.
+ ///
+ /// Sets the capacity map.
+ /// \return \c (*this)
+ Preflow& capacityMap(const CapacityMap& map) {
+ _capacity = ↦
+ return *this;
}
- /// Sets the edge map of the capacities to _cap.
+ /// \brief Sets the flow map.
+ ///
+ /// Sets the flow map.
+ /// \return \c (*this)
+ Preflow& flowMap(FlowMap& map) {
+ if (_local_flow) {
+ delete _flow;
+ _local_flow = false;
+ }
+ _flow = ↦
+ return *this;
+ }
- /// Sets the edge map of the capacities to _cap.
- ///
- void capacityMap(const CapacityMap& _cap) {
- _capacity=&_cap;
- status=AFTER_NOTHING;
+ /// \brief Returns the flow map.
+ ///
+ /// \return The flow map.
+ const FlowMap& flowMap() {
+ return *_flow;
}
- /// Returns a reference to capacity map.
- /// Returns a reference to capacity map.
- ///
- const CapacityMap &capacityMap() const {
- return *_capacity;
+ /// \brief Sets the elevator.
+ ///
+ /// Sets the elevator.
+ /// \return \c (*this)
+ Preflow& elevator(Elevator& elevator) {
+ if (_local_level) {
+ delete _level;
+ _local_level = false;
+ }
+ _level = &elevator;
+ return *this;
}
- /// Sets the edge map of the flows to _flow.
+ /// \brief Returns the elevator.
+ ///
+ /// \return The elevator.
+ const Elevator& elevator() {
+ return *_level;
+ }
- /// Sets the edge map of the flows to _flow.
- ///
- void flowMap(FlowMap& _f) {
- _flow=&_f;
- flow_prop=NO_FLOW;
- status=AFTER_NOTHING;
+ /// \brief Sets the source node.
+ ///
+ /// Sets the source node.
+ /// \return \c (*this)
+ Preflow& source(const Node& node) {
+ _source = node;
+ return *this;
}
-
- /// Returns a reference to flow map.
- /// Returns a reference to flow map.
- ///
- const FlowMap &flowMap() const {
- return *_flow;
+ /// \brief Sets the target node.
+ ///
+ /// Sets the target node.
+ /// \return \c (*this)
+ Preflow& target(const Node& node) {
+ _target = node;
+ return *this;
}
+
+ /// \brief Sets the tolerance used by algorithm.
+ ///
+ /// Sets the tolerance used by algorithm.
+ Preflow& tolerance(const Tolerance& tolerance) const {
+ _tolerance = tolerance;
+ return *this;
+ }
- private:
+ /// \brief Returns the tolerance used by algorithm.
+ ///
+ /// Returns the tolerance used by algorithm.
+ const Tolerance& tolerance() const {
+ return tolerance;
+ }
+
+ /// \name Execution control The simplest way to execute the
+ /// algorithm is to use one of the member functions called \c
+ /// run().
+ /// \n
+ /// If you need more control on initial solution or
+ /// execution then you have to call one \ref init() function and then
+ /// the startFirstPhase() and if you need the startSecondPhase().
+
+ ///@{
- int push(Node w, NNMap& next, VecNode& first) {
+ /// \brief Initializes the internal data structures.
+ ///
+ /// Initializes the internal data structures.
+ ///
+ void init() {
+ createStructures();
- int lev=level[w];
- Num exc=excess[w];
- int newlevel=_node_num; //bound on the next level of w
-
- for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
- if ( !_surely.positive((*_capacity)[e] - (*_flow)[e])) continue;
- Node v=_g->target(e);
-
- if( lev > level[v] ) { //Push is allowed now
-
- if ( excess[v] == 0 && v!=_target && v!=_source ) {
- next.set(v,first[level[v]]);
- first[level[v]]=v;
- }
+ _phase = true;
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ _excess->set(n, 0);
+ }
- Num cap=(*_capacity)[e];
- Num flo=(*_flow)[e];
- Num remcap=cap-flo;
-
- if ( ! _surely.less(remcap, exc) ) { //A nonsaturating push.
-
- _flow->set(e, flo+exc);
- excess.set(v, excess[v]+exc);
- exc=0;
- break;
-
- } else { //A saturating push.
- _flow->set(e, cap);
- excess.set(v, excess[v]+remcap);
- exc-=remcap;
- }
- } else if ( newlevel > level[v] ) newlevel = level[v];
- } //for out edges wv
+ for (EdgeIt e(_graph); e != INVALID; ++e) {
+ _flow->set(e, 0);
+ }
- if ( exc != 0 ) {
- for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
-
- if ( !_surely.positive((*_flow)[e]) ) continue;
- Node v=_g->source(e);
-
- if( lev > level[v] ) { //Push is allowed now
+ typename Graph::template NodeMap<bool> reached(_graph, false);
+
+ _level->initStart();
+ _level->initAddItem(_target);
- if ( excess[v] == 0 && v!=_target && v!=_source ) {
- next.set(v,first[level[v]]);
- first[level[v]]=v;
+ std::vector<Node> queue;
+ reached.set(_source, true);
+
+ queue.push_back(_target);
+ reached.set(_target, true);
+ while (!queue.empty()) {
+ std::vector<Node> nqueue;
+ for (int i = 0; i < int(queue.size()); ++i) {
+ Node n = queue[i];
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Node u = _graph.source(e);
+ if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
+ reached.set(u, true);
+ _level->initAddItem(u);
+ nqueue.push_back(u);
}
+ }
+ }
+ queue.swap(nqueue);
+ }
+ _level->initFinish();
- Num flo=(*_flow)[e];
+ for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
+ if (_tolerance.positive((*_capacity)[e])) {
+ Node u = _graph.target(e);
+ if ((*_level)[u] == _level->maxLevel()) continue;
+ _flow->set(e, (*_capacity)[e]);
+ _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
+ if (u != _target && !_level->active(u)) {
+ _level->activate(u);
+ }
+ }
+ }
+ }
- if ( !_surely.less(flo, exc) ) { //A nonsaturating push.
+ /// \brief Initializes the internal data structures.
+ ///
+ /// Initializes the internal data structures and sets the initial
+ /// flow to the given \c flowMap. The \c flowMap should contain a
+ /// flow or at least a preflow, ie. in each node excluding the
+ /// target the incoming flow should greater or equal to the
+ /// outgoing flow.
+ /// \return %False when the given \c flowMap is not a preflow.
+ template <typename FlowMap>
+ bool flowInit(const FlowMap& flowMap) {
+ createStructures();
- _flow->set(e, flo-exc);
- excess.set(v, excess[v]+exc);
- exc=0;
- break;
- } else { //A saturating push.
+ for (EdgeIt e(_graph); e != INVALID; ++e) {
+ _flow->set(e, flowMap[e]);
+ }
- excess.set(v, excess[v]+flo);
- exc-=flo;
- _flow->set(e,0);
- }
- } else if ( newlevel > level[v] ) newlevel = level[v];
- } //for in edges vw
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ Value excess = 0;
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+ excess += (*_flow)[e];
+ }
+ for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+ excess -= (*_flow)[e];
+ }
+ if (excess < 0 && n != _source) return false;
+ _excess->set(n, excess);
+ }
- } // if w still has excess after the out edge for cycle
+ typename Graph::template NodeMap<bool> reached(_graph, false);
- excess.set(w, exc);
-
- return newlevel;
- }
-
-
-
- void preflowPreproc(VecNode& first, NNMap& next,
- VecNode& level_list, NNMap& left, NNMap& right)
- {
- for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num);
- std::queue<Node> bfs_queue;
-
- if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
- //Reverse_bfs from t in the residual graph,
- //to find the starting level.
- level.set(_target,0);
- bfs_queue.push(_target);
-
- while ( !bfs_queue.empty() ) {
-
- Node v=bfs_queue.front();
- bfs_queue.pop();
- int l=level[v]+1;
-
- for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
- if ( !_surely.positive((*_capacity)[e] - (*_flow)[e] )) continue;
- Node w=_g->source(e);
- if ( level[w] == _node_num && w != _source ) {
- bfs_queue.push(w);
- Node z=level_list[l];
- if ( z!=INVALID ) left.set(z,w);
- right.set(w,z);
- level_list[l]=w;
- level.set(w, l);
+ _level->initStart();
+ _level->initAddItem(_target);
+
+ std::vector<Node> queue;
+ reached.set(_source, true);
+
+ queue.push_back(_target);
+ reached.set(_target, true);
+ while (!queue.empty()) {
+ _level->initNewLevel();
+ std::vector<Node> nqueue;
+ for (int i = 0; i < int(queue.size()); ++i) {
+ Node n = queue[i];
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Node u = _graph.source(e);
+ if (!reached[u] &&
+ _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
+ reached.set(u, true);
+ _level->initAddItem(u);
+ nqueue.push_back(u);
}
}
-
- for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
- if ( !_surely.positive((*_flow)[e]) ) continue;
- Node w=_g->target(e);
- if ( level[w] == _node_num && w != _source ) {
- bfs_queue.push(w);
- Node z=level_list[l];
- if ( z!=INVALID ) left.set(z,w);
- right.set(w,z);
- level_list[l]=w;
- level.set(w, l);
- }
- }
- } //while
- } //if
-
-
- switch (flow_prop) {
- case NO_FLOW:
- for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0);
- case ZERO_FLOW:
- for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
-
- //Reverse_bfs from t, to find the starting level.
- level.set(_target,0);
- bfs_queue.push(_target);
-
- while ( !bfs_queue.empty() ) {
-
- Node v=bfs_queue.front();
- bfs_queue.pop();
- int l=level[v]+1;
-
- for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
- Node w=_g->source(e);
- if ( level[w] == _node_num && w != _source ) {
- bfs_queue.push(w);
- Node z=level_list[l];
- if ( z!=INVALID ) left.set(z,w);
- right.set(w,z);
- level_list[l]=w;
- level.set(w, l);
+ for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Node v = _graph.target(e);
+ if (!reached[v] && _tolerance.positive((*_flow)[e])) {
+ reached.set(v, true);
+ _level->initAddItem(v);
+ nqueue.push_back(v);
}
}
}
-
- //the starting flow
- for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
- Num c=(*_capacity)[e];
- if ( !_surely.positive(c) ) continue;
- Node w=_g->target(e);
- if ( level[w] < _node_num ) {
- if ( excess[w] == 0 && w!=_target ) { //putting into the stack
- next.set(w,first[level[w]]);
- first[level[w]]=w;
- }
- _flow->set(e, c);
- excess.set(w, excess[w]+c);
- }
- }
- break;
-
- case GEN_FLOW:
- for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
- {
- Num exc=0;
- for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e];
- for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e];
- if (!_surely.positive(exc)) {
- exc = 0;
- }
- excess.set(_target,exc);
- }
-
- //the starting flow
- for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e) {
- Num rem=(*_capacity)[e]-(*_flow)[e];
- if ( !_surely.positive(rem) ) continue;
- Node w=_g->target(e);
- if ( level[w] < _node_num ) {
- if ( excess[w] == 0 && w!=_target ) { //putting into the stack
- next.set(w,first[level[w]]);
- first[level[w]]=w;
- }
- _flow->set(e, (*_capacity)[e]);
- excess.set(w, excess[w]+rem);
+ queue.swap(nqueue);
+ }
+ _level->initFinish();
+
+ for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
+ Value rem = (*_capacity)[e] - (*_flow)[e];
+ if (_tolerance.positive(rem)) {
+ Node u = _graph.target(e);
+ if ((*_level)[u] == _level->maxLevel()) continue;
+ _flow->set(e, (*_capacity)[e]);
+ _excess->set(u, (*_excess)[u] + rem);
+ if (u != _target && !_level->active(u)) {
+ _level->activate(u);
}
}
-
- for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) {
- if ( !_surely.positive((*_flow)[e]) ) continue;
- Node w=_g->source(e);
- if ( level[w] < _node_num ) {
- if ( excess[w] == 0 && w!=_target ) {
- next.set(w,first[level[w]]);
- first[level[w]]=w;
- }
- excess.set(w, excess[w]+(*_flow)[e]);
- _flow->set(e, 0);
- }
- }
- break;
-
- case PRE_FLOW:
- //the starting flow
- for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
- Num rem=(*_capacity)[e]-(*_flow)[e];
- if ( !_surely.positive(rem) ) continue;
- Node w=_g->target(e);
- if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]);
+ }
+ for (InEdgeIt e(_graph, _source); e != INVALID; ++e) {
+ Value rem = (*_flow)[e];
+ if (_tolerance.positive(rem)) {
+ Node v = _graph.source(e);
+ if ((*_level)[v] == _level->maxLevel()) continue;
+ _flow->set(e, 0);
+ _excess->set(v, (*_excess)[v] + rem);
+ if (v != _target && !_level->active(v)) {
+ _level->activate(v);
+ }
}
-
- for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
- if ( !_surely.positive((*_flow)[e]) ) continue;
- Node w=_g->source(e);
- if ( level[w] < _node_num ) _flow->set(e, 0);
+ }
+ return true;
+ }
+
+ /// \brief Starts the first phase of the preflow algorithm.
+ ///
+ /// The preflow algorithm consists of two phases, this method runs
+ /// the first phase. After the first phase the maximum flow value
+ /// and a minimum value cut can already be computed, although a
+ /// maximum flow is not yet obtained. So after calling this method
+ /// \ref flowValue() returns the value of a maximum flow and \ref
+ /// minCut() returns a minimum cut.
+ /// \pre One of the \ref init() functions should be called.
+ void startFirstPhase() {
+ _phase = true;
+
+ Node n = _level->highestActive();
+ int level = _level->highestActiveLevel();
+ while (n != INVALID) {
+ int num = _node_num;
+
+ while (num > 0 && n != INVALID) {
+ Value excess = (*_excess)[n];
+ int new_level = _level->maxLevel();
+
+ for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_capacity)[e] - (*_flow)[e];
+ if (!_tolerance.positive(rem)) continue;
+ Node v = _graph.target(e);
+ if ((*_level)[v] < level) {
+ if (!_level->active(v) && v != _target) {
+ _level->activate(v);
+ }
+ if (!_tolerance.less(rem, excess)) {
+ _flow->set(e, (*_flow)[e] + excess);
+ _excess->set(v, (*_excess)[v] + excess);
+ excess = 0;
+ goto no_more_push_1;
+ } else {
+ excess -= rem;
+ _excess->set(v, (*_excess)[v] + rem);
+ _flow->set(e, (*_capacity)[e]);
+ }
+ } else if (new_level > (*_level)[v]) {
+ new_level = (*_level)[v];
+ }
+ }
+
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_flow)[e];
+ if (!_tolerance.positive(rem)) continue;
+ Node v = _graph.source(e);
+ if ((*_level)[v] < level) {
+ if (!_level->active(v) && v != _target) {
+ _level->activate(v);
+ }
+ if (!_tolerance.less(rem, excess)) {
+ _flow->set(e, (*_flow)[e] - excess);
+ _excess->set(v, (*_excess)[v] + excess);
+ excess = 0;
+ goto no_more_push_1;
+ } else {
+ excess -= rem;
+ _excess->set(v, (*_excess)[v] + rem);
+ _flow->set(e, 0);
+ }
+ } else if (new_level > (*_level)[v]) {
+ new_level = (*_level)[v];
+ }
+ }
+
+ no_more_push_1:
+
+ _excess->set(n, excess);
+
+ if (excess != 0) {
+ if (new_level + 1 < _level->maxLevel()) {
+ _level->liftHighestActive(new_level + 1);
+ } else {
+ _level->liftHighestActiveToTop();
+ }
+ if (_level->emptyLevel(level)) {
+ _level->liftToTop(level);
+ }
+ } else {
+ _level->deactivate(n);
+ }
+
+ n = _level->highestActive();
+ level = _level->highestActiveLevel();
+ --num;
}
- //computing the excess
- for(NodeIt w(*_g); w!=INVALID; ++w) {
- Num exc=0;
- for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e];
- for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e];
- if (!_surely.positive(exc)) {
- exc = 0;
- }
- excess.set(w,exc);
-
- //putting the active nodes into the stack
- int lev=level[w];
- if ( exc != 0 && lev < _node_num && Node(w) != _target ) {
- next.set(w,first[lev]);
- first[lev]=w;
+ num = _node_num * 20;
+ while (num > 0 && n != INVALID) {
+ Value excess = (*_excess)[n];
+ int new_level = _level->maxLevel();
+
+ for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_capacity)[e] - (*_flow)[e];
+ if (!_tolerance.positive(rem)) continue;
+ Node v = _graph.target(e);
+ if ((*_level)[v] < level) {
+ if (!_level->active(v) && v != _target) {
+ _level->activate(v);
+ }
+ if (!_tolerance.less(rem, excess)) {
+ _flow->set(e, (*_flow)[e] + excess);
+ _excess->set(v, (*_excess)[v] + excess);
+ excess = 0;
+ goto no_more_push_2;
+ } else {
+ excess -= rem;
+ _excess->set(v, (*_excess)[v] + rem);
+ _flow->set(e, (*_capacity)[e]);
+ }
+ } else if (new_level > (*_level)[v]) {
+ new_level = (*_level)[v];
}
- }
- break;
- } //switch
- } //preflowPreproc
+ }
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_flow)[e];
+ if (!_tolerance.positive(rem)) continue;
+ Node v = _graph.source(e);
+ if ((*_level)[v] < level) {
+ if (!_level->active(v) && v != _target) {
+ _level->activate(v);
+ }
+ if (!_tolerance.less(rem, excess)) {
+ _flow->set(e, (*_flow)[e] - excess);
+ _excess->set(v, (*_excess)[v] + excess);
+ excess = 0;
+ goto no_more_push_2;
+ } else {
+ excess -= rem;
+ _excess->set(v, (*_excess)[v] + rem);
+ _flow->set(e, 0);
+ }
+ } else if (new_level > (*_level)[v]) {
+ new_level = (*_level)[v];
+ }
+ }
- void relabel(Node w, int newlevel, VecNode& first, NNMap& next,
- VecNode& level_list, NNMap& left,
- NNMap& right, int& b, int& k, bool what_heur )
- {
+ no_more_push_2:
- int lev=level[w];
+ _excess->set(n, excess);
- Node right_n=right[w];
- Node left_n=left[w];
+ if (excess != 0) {
+ if (new_level + 1 < _level->maxLevel()) {
+ _level->liftActiveOn(level, new_level + 1);
+ } else {
+ _level->liftActiveToTop(level);
+ }
+ if (_level->emptyLevel(level)) {
+ _level->liftToTop(level);
+ }
+ } else {
+ _level->deactivate(n);
+ }
- //unlacing starts
- if ( right_n!=INVALID ) {
- if ( left_n!=INVALID ) {
- right.set(left_n, right_n);
- left.set(right_n, left_n);
- } else {
- level_list[lev]=right_n;
- left.set(right_n, INVALID);
+ while (level >= 0 && _level->activeFree(level)) {
+ --level;
+ }
+ if (level == -1) {
+ n = _level->highestActive();
+ level = _level->highestActiveLevel();
+ } else {
+ n = _level->activeOn(level);
+ }
+ --num;
}
- } else {
- if ( left_n!=INVALID ) {
- right.set(left_n, INVALID);
- } else {
- level_list[lev]=INVALID;
+ }
+ }
+
+ /// \brief Starts the second phase of the preflow algorithm.
+ ///
+ /// The preflow algorithm consists of two phases, this method runs
+ /// the second phase. After calling \ref init() and \ref
+ /// startFirstPhase() and then \ref startSecondPhase(), \ref
+ /// flowMap() return a maximum flow, \ref flowValue() returns the
+ /// value of a maximum flow, \ref minCut() returns a minimum cut
+ /// \pre The \ref init() and startFirstPhase() functions should be
+ /// called before.
+ void startSecondPhase() {
+ _phase = false;
+
+ typename Graph::template NodeMap<bool> reached(_graph, false);
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ reached.set(n, (*_level)[n] < _level->maxLevel());
+ }
+
+ _level->initStart();
+ _level->initAddItem(_source);
+
+ std::vector<Node> queue;
+ queue.push_back(_source);
+ reached.set(_source, true);
+
+ while (!queue.empty()) {
+ _level->initNewLevel();
+ std::vector<Node> nqueue;
+ for (int i = 0; i < int(queue.size()); ++i) {
+ Node n = queue[i];
+ for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Node v = _graph.target(e);
+ if (!reached[v] && _tolerance.positive((*_flow)[e])) {
+ reached.set(v, true);
+ _level->initAddItem(v);
+ nqueue.push_back(v);
+ }
+ }
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Node u = _graph.source(e);
+ if (!reached[u] &&
+ _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
+ reached.set(u, true);
+ _level->initAddItem(u);
+ nqueue.push_back(u);
+ }
+ }
}
+ queue.swap(nqueue);
}
- //unlacing ends
+ _level->initFinish();
- if ( level_list[lev]==INVALID ) {
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ if ((*_excess)[n] > 0 && _target != n) {
+ _level->activate(n);
+ }
+ }
- //gapping starts
- for (int i=lev; i!=k ; ) {
- Node v=level_list[++i];
- while ( v!=INVALID ) {
- level.set(v,_node_num);
- v=right[v];
+ Node n;
+ while ((n = _level->highestActive()) != INVALID) {
+ Value excess = (*_excess)[n];
+ int level = _level->highestActiveLevel();
+ int new_level = _level->maxLevel();
+
+ for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_capacity)[e] - (*_flow)[e];
+ if (!_tolerance.positive(rem)) continue;
+ Node v = _graph.target(e);
+ if ((*_level)[v] < level) {
+ if (!_level->active(v) && v != _source) {
+ _level->activate(v);
+ }
+ if (!_tolerance.less(rem, excess)) {
+ _flow->set(e, (*_flow)[e] + excess);
+ _excess->set(v, (*_excess)[v] + excess);
+ excess = 0;
+ goto no_more_push;
+ } else {
+ excess -= rem;
+ _excess->set(v, (*_excess)[v] + rem);
+ _flow->set(e, (*_capacity)[e]);
+ }
+ } else if (new_level > (*_level)[v]) {
+ new_level = (*_level)[v];
}
- level_list[i]=INVALID;
- if ( !what_heur ) first[i]=INVALID;
}
- level.set(w,_node_num);
- b=lev-1;
- k=b;
- //gapping ends
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_flow)[e];
+ if (!_tolerance.positive(rem)) continue;
+ Node v = _graph.source(e);
+ if ((*_level)[v] < level) {
+ if (!_level->active(v) && v != _source) {
+ _level->activate(v);
+ }
+ if (!_tolerance.less(rem, excess)) {
+ _flow->set(e, (*_flow)[e] - excess);
+ _excess->set(v, (*_excess)[v] + excess);
+ excess = 0;
+ goto no_more_push;
+ } else {
+ excess -= rem;
+ _excess->set(v, (*_excess)[v] + rem);
+ _flow->set(e, 0);
+ }
+ } else if (new_level > (*_level)[v]) {
+ new_level = (*_level)[v];
+ }
+ }
- } else {
+ no_more_push:
- if ( newlevel == _node_num ) level.set(w,_node_num);
- else {
- level.set(w,++newlevel);
- next.set(w,first[newlevel]);
- first[newlevel]=w;
- if ( what_heur ) b=newlevel;
- if ( k < newlevel ) ++k; //now k=newlevel
- Node z=level_list[newlevel];
- if ( z!=INVALID ) left.set(z,w);
- right.set(w,z);
- left.set(w,INVALID);
- level_list[newlevel]=w;
+ _excess->set(n, excess);
+
+ if (excess != 0) {
+ if (new_level + 1 < _level->maxLevel()) {
+ _level->liftHighestActive(new_level + 1);
+ } else {
+ // Calculation error
+ _level->liftHighestActiveToTop();
+ }
+ if (_level->emptyLevel(level)) {
+ // Calculation error
+ _level->liftToTop(level);
+ }
+ } else {
+ _level->deactivate(n);
}
+
}
- } //relabel
+ }
- };
+ /// \brief Runs the preflow algorithm.
+ ///
+ /// Runs the preflow algorithm.
+ /// \note pf.run() is just a shortcut of the following code.
+ /// \code
+ /// pf.init();
+ /// pf.startFirstPhase();
+ /// pf.startSecondPhase();
+ /// \endcode
+ void run() {
+ init();
+ startFirstPhase();
+ startSecondPhase();
+ }
- ///\ingroup max_flow
- ///\brief Function type interface for Preflow algorithm.
- ///
- ///Function type interface for Preflow algorithm.
- ///\sa Preflow
- template<class GR, class CM, class FM>
- Preflow<GR,typename CM::Value,CM,FM> preflow(const GR &g,
- typename GR::Node source,
- typename GR::Node target,
- const CM &cap,
- FM &flow
- )
- {
- return Preflow<GR,typename CM::Value,CM,FM>(g,source,target,cap,flow);
- }
+ /// \brief Runs the preflow algorithm to compute the minimum cut.
+ ///
+ /// Runs the preflow algorithm to compute the minimum cut.
+ /// \note pf.runMinCut() is just a shortcut of the following code.
+ /// \code
+ /// pf.init();
+ /// pf.startFirstPhase();
+ /// \endcode
+ void runMinCut() {
+ init();
+ startFirstPhase();
+ }
+
+ /// @}
+
+ /// \name Query Functions
+ /// The result of the %Dijkstra algorithm can be obtained using these
+ /// functions.\n
+ /// Before the use of these functions,
+ /// either run() or start() must be called.
+
+ ///@{
+
+ /// \brief Returns the value of the maximum flow.
+ ///
+ /// Returns the value of the maximum flow by returning the excess
+ /// of the target node \c t. This value equals to the value of
+ /// the maximum flow already after the first phase.
+ Value flowValue() const {
+ return (*_excess)[_target];
+ }
+
+ /// \brief Returns true when the node is on the source side of minimum cut.
+ ///
+ /// Returns true when the node is on the source side of minimum
+ /// cut. This method can be called both after running \ref
+ /// startFirstPhase() and \ref startSecondPhase().
+ bool minCut(const Node& node) const {
+ return ((*_level)[node] == _level->maxLevel()) == _phase;
+ }
+
+ /// \brief Returns a minimum value cut.
+ ///
+ /// Sets the \c cutMap to the characteristic vector of a minimum value
+ /// cut. This method can be called both after running \ref
+ /// startFirstPhase() and \ref startSecondPhase(). The result after second
+ /// phase could be changed slightly if inexact computation is used.
+ /// \pre The \c cutMap should be a bool-valued node-map.
+ template <typename CutMap>
+ void minCutMap(CutMap& cutMap) const {
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ cutMap.set(n, minCut(n));
+ }
+ }
-} //namespace lemon
+ /// \brief Returns the flow on the edge.
+ ///
+ /// Sets the \c flowMap to the flow on the edges. This method can
+ /// be called after the second phase of algorithm.
+ Value flow(const Edge& edge) const {
+ return (*_flow)[edge];
+ }
+
+ /// @}
+ };
+}
-#endif //LEMON_PREFLOW_H
+#endif
Modified: lemon/trunk/test/preflow_test.cc
==============================================================================
--- lemon/trunk/test/preflow_test.cc (original)
+++ lemon/trunk/test/preflow_test.cc Sat Nov 17 21:58:11 2007
@@ -28,7 +28,7 @@
using namespace lemon;
-void check_Preflow()
+void checkPreflow()
{
typedef int VType;
typedef concepts::Graph Graph;
@@ -37,35 +37,40 @@
typedef Graph::Edge Edge;
typedef concepts::ReadMap<Edge,VType> CapMap;
typedef concepts::ReadWriteMap<Edge,VType> FlowMap;
- typedef concepts::ReadWriteMap<Node,bool> CutMap;
+ typedef concepts::WriteMap<Node,bool> CutMap;
- typedef Preflow<Graph, int, CapMap, FlowMap> PType;
-
Graph g;
Node n;
+ Edge e;
CapMap cap;
FlowMap flow;
CutMap cut;
- PType preflow_test(g,n,n,cap,flow);
+ Preflow<Graph, CapMap>::DefFlowMap<FlowMap>::Create preflow_test(g,cap,n,n);
- preflow_test.run();
- preflow_test.flowValue();
- preflow_test.source(n);
+ preflow_test.capacityMap(cap);
+ flow = preflow_test.flowMap();
preflow_test.flowMap(flow);
+ preflow_test.source(n);
+ preflow_test.target(n);
+
+ preflow_test.init();
+ preflow_test.flowInit(cap);
+ preflow_test.startFirstPhase();
+ preflow_test.startSecondPhase();
+ preflow_test.run();
+ preflow_test.runMinCut();
- preflow_test.phase1(PType::NO_FLOW);
- preflow_test.minCut(cut);
+ preflow_test.flowValue();
+ preflow_test.minCut(n);
+ preflow_test.minCutMap(cut);
+ preflow_test.flow(e);
- preflow_test.phase2();
- preflow_test.target(n);
- preflow_test.capacityMap(cap);
- preflow_test.minMinCut(cut);
- preflow_test.maxMinCut(cut);
}
-int cut_value ( SmartGraph& g, SmartGraph::NodeMap<bool>& cut,
- SmartGraph::EdgeMap<int>& cap) {
+int cutValue (const SmartGraph& g,
+ const SmartGraph::NodeMap<bool>& cut,
+ const SmartGraph::EdgeMap<int>& cap) {
int c=0;
for(SmartGraph::EdgeIt e(g); e!=INVALID; ++e) {
@@ -74,6 +79,29 @@
return c;
}
+bool checkFlow(const SmartGraph& g,
+ const SmartGraph::EdgeMap<int>& flow,
+ const SmartGraph::EdgeMap<int>& cap,
+ SmartGraph::Node s, SmartGraph::Node t) {
+
+ for (SmartGraph::EdgeIt e(g); e != INVALID; ++e) {
+ if (flow[e] < 0 || flow[e] > cap[e]) return false;
+ }
+
+ for (SmartGraph::NodeIt n(g); n != INVALID; ++n) {
+ if (n == s || n == t) continue;
+ int sum = 0;
+ for (SmartGraph::OutEdgeIt e(g, n); e != INVALID; ++e) {
+ sum += flow[e];
+ }
+ for (SmartGraph::InEdgeIt e(g, n); e != INVALID; ++e) {
+ sum -= flow[e];
+ }
+ if (sum != 0) return false;
+ }
+ return true;
+}
+
int main() {
typedef SmartGraph Graph;
@@ -85,7 +113,7 @@
typedef Graph::EdgeMap<int> FlowMap;
typedef Graph::NodeMap<bool> CutMap;
- typedef Preflow<Graph, int> PType;
+ typedef Preflow<Graph, CapMap> PType;
std::string f_name;
if( getenv("srcdir") )
@@ -102,74 +130,50 @@
CapMap cap(g);
readDimacs(file, g, cap, s, t);
- FlowMap flow(g,0);
-
-
-
- PType preflow_test(g, s, t, cap, flow);
- preflow_test.run(PType::ZERO_FLOW);
+ PType preflow_test(g, cap, s, t);
+ preflow_test.run();
- CutMap min_cut(g,false);
- preflow_test.minCut(min_cut);
- int min_cut_value=cut_value(g,min_cut,cap);
-
- CutMap min_min_cut(g,false);
- preflow_test.minMinCut(min_min_cut);
- int min_min_cut_value=cut_value(g,min_min_cut,cap);
-
- CutMap max_min_cut(g,false);
- preflow_test.maxMinCut(max_min_cut);
- int max_min_cut_value=cut_value(g,max_min_cut,cap);
+ check(checkFlow(g, preflow_test.flowMap(), cap, s, t),
+ "The flow is not feasible.");
- check(preflow_test.flowValue() == min_cut_value &&
- min_cut_value == min_min_cut_value &&
- min_min_cut_value == max_min_cut_value,
+ CutMap min_cut(g);
+ preflow_test.minCutMap(min_cut);
+ int min_cut_value=cutValue(g,min_cut,cap);
+
+ check(preflow_test.flowValue() == min_cut_value,
"The max flow value is not equal to the three min cut values.");
- int flow_value=preflow_test.flowValue();
-
+ FlowMap flow(g);
+ flow = preflow_test.flowMap();
+ int flow_value=preflow_test.flowValue();
for(EdgeIt e(g); e!=INVALID; ++e) cap[e]=2*cap[e];
- preflow_test.capacityMap(cap);
-
- preflow_test.phase1(PType::PRE_FLOW);
+ preflow_test.flowInit(flow);
+ preflow_test.startFirstPhase();
- CutMap min_cut1(g,false);
- preflow_test.minCut(min_cut1);
- min_cut_value=cut_value(g,min_cut1,cap);
+ CutMap min_cut1(g);
+ preflow_test.minCutMap(min_cut1);
+ min_cut_value=cutValue(g,min_cut1,cap);
check(preflow_test.flowValue() == min_cut_value &&
min_cut_value == 2*flow_value,
"The max flow value or the min cut value is wrong.");
- preflow_test.phase2();
+ preflow_test.startSecondPhase();
- CutMap min_cut2(g,false);
- preflow_test.minCut(min_cut2);
- min_cut_value=cut_value(g,min_cut2,cap);
-
- CutMap min_min_cut2(g,false);
- preflow_test.minMinCut(min_min_cut2);
- min_min_cut_value=cut_value(g,min_min_cut2,cap);
-
- preflow_test.maxMinCut(max_min_cut);
- max_min_cut_value=cut_value(g,max_min_cut,cap);
+ check(checkFlow(g, preflow_test.flowMap(), cap, s, t),
+ "The flow is not feasible.");
+ CutMap min_cut2(g);
+ preflow_test.minCutMap(min_cut2);
+ min_cut_value=cutValue(g,min_cut2,cap);
+
check(preflow_test.flowValue() == min_cut_value &&
- min_cut_value == min_min_cut_value &&
- min_min_cut_value == max_min_cut_value &&
min_cut_value == 2*flow_value,
"The max flow value or the three min cut values were not doubled");
-
- EdgeIt e(g);
- for( int i=1; i==10; ++i ) {
- flow.set(e,0);
- ++e;
- }
-
preflow_test.flowMap(flow);
NodeIt tmp1(g,s);
@@ -185,19 +189,13 @@
preflow_test.run();
- CutMap min_cut3(g,false);
- preflow_test.minCut(min_cut3);
- min_cut_value=cut_value(g,min_cut3,cap);
-
- CutMap min_min_cut3(g,false);
- preflow_test.minMinCut(min_min_cut3);
- min_min_cut_value=cut_value(g,min_min_cut3,cap);
+ CutMap min_cut3(g);
+ preflow_test.minCutMap(min_cut3);
+ min_cut_value=cutValue(g,min_cut3,cap);
- preflow_test.maxMinCut(max_min_cut);
- max_min_cut_value=cut_value(g,max_min_cut,cap);
- check(preflow_test.flowValue() == min_cut_value &&
- min_cut_value == min_min_cut_value &&
- min_min_cut_value == max_min_cut_value,
+ check(preflow_test.flowValue() == min_cut_value,
"The max flow value or the three min cut values are incorrect.");
+
+ return 0;
}
More information about the Lemon-commits
mailing list