lemon: lemon/preflow.h@2b5496c62ccd (annotated)

alpar@404	1	/* -- mode: C++; indent-tabs-mode: nil; --
alpar@404	2	*
alpar@404	3	* This file is a part of LEMON, a generic C++ optimization library.
alpar@404	4	*
alpar@463	5	* Copyright (C) 2003-2009
alpar@404	6	* Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
alpar@404	7	* (Egervary Research Group on Combinatorial Optimization, EGRES).
alpar@404	8	*
alpar@404	9	* Permission to use, modify and distribute this software is granted
alpar@404	10	* provided that this copyright notice appears in all copies. For
alpar@404	11	* precise terms see the accompanying LICENSE file.
alpar@404	12	*
alpar@404	13	* This software is provided "AS IS" with no warranty of any kind,
alpar@404	14	* express or implied, and with no claim as to its suitability for any
alpar@404	15	* purpose.
alpar@404	16	*
alpar@404	17	*/
alpar@404	18
alpar@404	19	#ifndef LEMON_PREFLOW_H
alpar@404	20	#define LEMON_PREFLOW_H
alpar@404	21
alpar@404	22	#include <lemon/tolerance.h>
alpar@404	23	#include <lemon/elevator.h>
alpar@404	24
alpar@404	25	/// \file
alpar@404	26	/// \ingroup max_flow
alpar@404	27	/// \brief Implementation of the preflow algorithm.
alpar@404	28
alpar@404	29	namespace lemon {
alpar@404	30
alpar@404	31	/// \brief Default traits class of Preflow class.
alpar@404	32	///
alpar@404	33	/// Default traits class of Preflow class.
kpeter@408	34	/// \tparam _Digraph Digraph type.
kpeter@408	35	/// \tparam _CapacityMap Capacity map type.
kpeter@408	36	template <typename _Digraph, typename _CapacityMap>
alpar@404	37	struct PreflowDefaultTraits {
alpar@404	38
kpeter@408	39	/// \brief The type of the digraph the algorithm runs on.
kpeter@408	40	typedef _Digraph Digraph;
alpar@404	41
alpar@404	42	/// \brief The type of the map that stores the arc capacities.
alpar@404	43	///
alpar@404	44	/// The type of the map that stores the arc capacities.
alpar@404	45	/// It must meet the \ref concepts::ReadMap "ReadMap" concept.
alpar@404	46	typedef _CapacityMap CapacityMap;
alpar@404	47
kpeter@408	48	/// \brief The type of the flow values.
alpar@404	49	typedef typename CapacityMap::Value Value;
alpar@404	50
kpeter@408	51	/// \brief The type of the map that stores the flow values.
alpar@404	52	///
kpeter@408	53	/// The type of the map that stores the flow values.
alpar@404	54	/// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
alpar@404	55	typedef typename Digraph::template ArcMap<Value> FlowMap;
alpar@404	56
alpar@404	57	/// \brief Instantiates a FlowMap.
alpar@404	58	///
alpar@404	59	/// This function instantiates a \ref FlowMap.
alpar@404	60	/// \param digraph The digraph, to which we would like to define
alpar@404	61	/// the flow map.
alpar@404	62	static FlowMap* createFlowMap(const Digraph& digraph) {
alpar@404	63	return new FlowMap(digraph);
alpar@404	64	}
alpar@404	65
kpeter@408	66	/// \brief The elevator type used by Preflow algorithm.
alpar@404	67	///
alpar@404	68	/// The elevator type used by Preflow algorithm.
alpar@404	69	///
alpar@404	70	/// \sa Elevator
alpar@404	71	/// \sa LinkedElevator
alpar@404	72	typedef LinkedElevator<Digraph, typename Digraph::Node> Elevator;
alpar@404	73
alpar@404	74	/// \brief Instantiates an Elevator.
alpar@404	75	///
kpeter@408	76	/// This function instantiates an \ref Elevator.
alpar@404	77	/// \param digraph The digraph, to which we would like to define
alpar@404	78	/// the elevator.
alpar@404	79	/// \param max_level The maximum level of the elevator.
alpar@404	80	static Elevator* createElevator(const Digraph& digraph, int max_level) {
alpar@404	81	return new Elevator(digraph, max_level);
alpar@404	82	}
alpar@404	83
alpar@404	84	/// \brief The tolerance used by the algorithm
alpar@404	85	///
alpar@404	86	/// The tolerance used by the algorithm to handle inexact computation.
alpar@404	87	typedef lemon::Tolerance<Value> Tolerance;
alpar@404	88
alpar@404	89	};
alpar@404	90
alpar@404	91
alpar@404	92	/// \ingroup max_flow
alpar@404	93	///
kpeter@408	94	/// \brief %Preflow algorithm class.
alpar@404	95	///
kpeter@408	96	/// This class provides an implementation of Goldberg-Tarjan's \e preflow
kpeter@408	97	/// \e push-relabel algorithm producing a flow of maximum value in a
kpeter@408	98	/// digraph. The preflow algorithms are the fastest known maximum
alpar@404	99	/// flow algorithms. The current implementation use a mixture of the
alpar@404	100	/// \e "highest label" and the \e "bound decrease" heuristics.
alpar@404	101	/// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$.
alpar@404	102	///
kpeter@408	103	/// The algorithm consists of two phases. After the first phase
kpeter@408	104	/// the maximum flow value and the minimum cut is obtained. The
kpeter@408	105	/// second phase constructs a feasible maximum flow on each arc.
alpar@404	106	///
kpeter@408	107	/// \tparam _Digraph The type of the digraph the algorithm runs on.
kpeter@408	108	/// \tparam _CapacityMap The type of the capacity map. The default map
kpeter@408	109	/// type is \ref concepts::Digraph::ArcMap "_Digraph::ArcMap<int>".
alpar@404	110	#ifdef DOXYGEN
kpeter@408	111	template <typename _Digraph, typename _CapacityMap, typename _Traits>
alpar@404	112	#else
kpeter@408	113	template <typename _Digraph,
kpeter@408	114	typename _CapacityMap = typename _Digraph::template ArcMap<int>,
kpeter@408	115	typename _Traits = PreflowDefaultTraits<_Digraph, _CapacityMap> >
alpar@404	116	#endif
alpar@404	117	class Preflow {
alpar@404	118	public:
alpar@404	119
kpeter@408	120	///The \ref PreflowDefaultTraits "traits class" of the algorithm.
alpar@404	121	typedef _Traits Traits;
kpeter@408	122	///The type of the digraph the algorithm runs on.
alpar@404	123	typedef typename Traits::Digraph Digraph;
kpeter@408	124	///The type of the capacity map.
alpar@404	125	typedef typename Traits::CapacityMap CapacityMap;
kpeter@408	126	///The type of the flow values.
alpar@404	127	typedef typename Traits::Value Value;
alpar@404	128
kpeter@408	129	///The type of the flow map.
alpar@404	130	typedef typename Traits::FlowMap FlowMap;
kpeter@408	131	///The type of the elevator.
alpar@404	132	typedef typename Traits::Elevator Elevator;
kpeter@408	133	///The type of the tolerance.
alpar@404	134	typedef typename Traits::Tolerance Tolerance;
alpar@404	135
alpar@404	136	private:
alpar@404	137
alpar@404	138	TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
alpar@404	139
alpar@404	140	const Digraph& _graph;
alpar@404	141	const CapacityMap* _capacity;
alpar@404	142
alpar@404	143	int _node_num;
alpar@404	144
alpar@404	145	Node _source, _target;
alpar@404	146
alpar@404	147	FlowMap* _flow;
alpar@404	148	bool _local_flow;
alpar@404	149
alpar@404	150	Elevator* _level;
alpar@404	151	bool _local_level;
alpar@404	152
alpar@404	153	typedef typename Digraph::template NodeMap<Value> ExcessMap;
alpar@404	154	ExcessMap* _excess;
alpar@404	155
alpar@404	156	Tolerance _tolerance;
alpar@404	157
alpar@404	158	bool _phase;
alpar@404	159
alpar@404	160
alpar@404	161	void createStructures() {
alpar@404	162	_node_num = countNodes(_graph);
alpar@404	163
alpar@404	164	if (!_flow) {
alpar@404	165	_flow = Traits::createFlowMap(_graph);
alpar@404	166	_local_flow = true;
alpar@404	167	}
alpar@404	168	if (!_level) {
alpar@404	169	_level = Traits::createElevator(_graph, _node_num);
alpar@404	170	_local_level = true;
alpar@404	171	}
alpar@404	172	if (!_excess) {
alpar@404	173	_excess = new ExcessMap(_graph);
alpar@404	174	}
alpar@404	175	}
alpar@404	176
alpar@404	177	void destroyStructures() {
alpar@404	178	if (_local_flow) {
alpar@404	179	delete _flow;
alpar@404	180	}
alpar@404	181	if (_local_level) {
alpar@404	182	delete _level;
alpar@404	183	}
alpar@404	184	if (_excess) {
alpar@404	185	delete _excess;
alpar@404	186	}
alpar@404	187	}
alpar@404	188
alpar@404	189	public:
alpar@404	190
alpar@404	191	typedef Preflow Create;
alpar@404	192
kpeter@408	193	///\name Named Template Parameters
alpar@404	194
alpar@404	195	///@{
alpar@404	196
alpar@404	197	template <typename _FlowMap>
alpar@406	198	struct SetFlowMapTraits : public Traits {
alpar@404	199	typedef _FlowMap FlowMap;
alpar@404	200	static FlowMap *createFlowMap(const Digraph&) {
alpar@405	201	LEMON_ASSERT(false, "FlowMap is not initialized");
alpar@405	202	return 0; // ignore warnings
alpar@404	203	}
alpar@404	204	};
alpar@404	205
alpar@404	206	/// \brief \ref named-templ-param "Named parameter" for setting
alpar@404	207	/// FlowMap type
alpar@404	208	///
alpar@404	209	/// \ref named-templ-param "Named parameter" for setting FlowMap
kpeter@408	210	/// type.
alpar@404	211	template <typename _FlowMap>
alpar@406	212	struct SetFlowMap
alpar@406	213	: public Preflow<Digraph, CapacityMap, SetFlowMapTraits<_FlowMap> > {
alpar@404	214	typedef Preflow<Digraph, CapacityMap,
alpar@406	215	SetFlowMapTraits<_FlowMap> > Create;
alpar@404	216	};
alpar@404	217
alpar@404	218	template <typename _Elevator>
alpar@406	219	struct SetElevatorTraits : public Traits {
alpar@404	220	typedef _Elevator Elevator;
alpar@404	221	static Elevator *createElevator(const Digraph&, int) {
alpar@405	222	LEMON_ASSERT(false, "Elevator is not initialized");
alpar@405	223	return 0; // ignore warnings
alpar@404	224	}
alpar@404	225	};
alpar@404	226
alpar@404	227	/// \brief \ref named-templ-param "Named parameter" for setting
alpar@404	228	/// Elevator type
alpar@404	229	///
alpar@404	230	/// \ref named-templ-param "Named parameter" for setting Elevator
kpeter@408	231	/// type. If this named parameter is used, then an external
kpeter@408	232	/// elevator object must be passed to the algorithm using the
kpeter@408	233	/// \ref elevator(Elevator&) "elevator()" function before calling
kpeter@408	234	/// \ref run() or \ref init().
kpeter@408	235	/// \sa SetStandardElevator
alpar@404	236	template <typename _Elevator>
alpar@406	237	struct SetElevator
alpar@406	238	: public Preflow<Digraph, CapacityMap, SetElevatorTraits<_Elevator> > {
alpar@404	239	typedef Preflow<Digraph, CapacityMap,
alpar@406	240	SetElevatorTraits<_Elevator> > Create;
alpar@404	241	};
alpar@404	242
alpar@404	243	template <typename _Elevator>
alpar@406	244	struct SetStandardElevatorTraits : public Traits {
alpar@404	245	typedef _Elevator Elevator;
alpar@404	246	static Elevator *createElevator(const Digraph& digraph, int max_level) {
alpar@404	247	return new Elevator(digraph, max_level);
alpar@404	248	}
alpar@404	249	};
alpar@404	250
alpar@404	251	/// \brief \ref named-templ-param "Named parameter" for setting
kpeter@408	252	/// Elevator type with automatic allocation
alpar@404	253	///
alpar@404	254	/// \ref named-templ-param "Named parameter" for setting Elevator
kpeter@408	255	/// type with automatic allocation.
kpeter@408	256	/// The Elevator should have standard constructor interface to be
kpeter@408	257	/// able to automatically created by the algorithm (i.e. the
kpeter@408	258	/// digraph and the maximum level should be passed to it).
kpeter@408	259	/// However an external elevator object could also be passed to the
kpeter@408	260	/// algorithm with the \ref elevator(Elevator&) "elevator()" function
kpeter@408	261	/// before calling \ref run() or \ref init().
kpeter@408	262	/// \sa SetElevator
alpar@404	263	template <typename _Elevator>
alpar@406	264	struct SetStandardElevator
alpar@404	265	: public Preflow<Digraph, CapacityMap,
alpar@406	266	SetStandardElevatorTraits<_Elevator> > {
alpar@404	267	typedef Preflow<Digraph, CapacityMap,
alpar@406	268	SetStandardElevatorTraits<_Elevator> > Create;
alpar@404	269	};
alpar@404	270
alpar@404	271	/// @}
alpar@404	272
alpar@404	273	protected:
alpar@404	274
alpar@404	275	Preflow() {}
alpar@404	276
alpar@404	277	public:
alpar@404	278
alpar@404	279
alpar@404	280	/// \brief The constructor of the class.
alpar@404	281	///
alpar@404	282	/// The constructor of the class.
alpar@404	283	/// \param digraph The digraph the algorithm runs on.
alpar@404	284	/// \param capacity The capacity of the arcs.
alpar@404	285	/// \param source The source node.
alpar@404	286	/// \param target The target node.
alpar@404	287	Preflow(const Digraph& digraph, const CapacityMap& capacity,
kpeter@408	288	Node source, Node target)
alpar@404	289	: _graph(digraph), _capacity(&capacity),
alpar@404	290	_node_num(0), _source(source), _target(target),
alpar@404	291	_flow(0), _local_flow(false),
alpar@404	292	_level(0), _local_level(false),
alpar@404	293	_excess(0), _tolerance(), _phase() {}
alpar@404	294
kpeter@408	295	/// \brief Destructor.
alpar@404	296	///
alpar@404	297	/// Destructor.
alpar@404	298	~Preflow() {
alpar@404	299	destroyStructures();
alpar@404	300	}
alpar@404	301
alpar@404	302	/// \brief Sets the capacity map.
alpar@404	303	///
alpar@404	304	/// Sets the capacity map.
kpeter@408	305	/// \return <tt>(*this)</tt>
alpar@404	306	Preflow& capacityMap(const CapacityMap& map) {
alpar@404	307	_capacity = &map;
alpar@404	308	return *this;
alpar@404	309	}
alpar@404	310
alpar@404	311	/// \brief Sets the flow map.
alpar@404	312	///
alpar@404	313	/// Sets the flow map.
kpeter@408	314	/// If you don't use this function before calling \ref run() or
kpeter@408	315	/// \ref init(), an instance will be allocated automatically.
kpeter@408	316	/// The destructor deallocates this automatically allocated map,
kpeter@408	317	/// of course.
kpeter@408	318	/// \return <tt>(*this)</tt>
alpar@404	319	Preflow& flowMap(FlowMap& map) {
alpar@404	320	if (_local_flow) {
alpar@404	321	delete _flow;
alpar@404	322	_local_flow = false;
alpar@404	323	}
alpar@404	324	_flow = &map;
alpar@404	325	return *this;
alpar@404	326	}
alpar@404	327
kpeter@408	328	/// \brief Sets the source node.
alpar@404	329	///
kpeter@408	330	/// Sets the source node.
kpeter@408	331	/// \return <tt>(*this)</tt>
kpeter@408	332	Preflow& source(const Node& node) {
kpeter@408	333	_source = node;
kpeter@408	334	return *this;
alpar@404	335	}
alpar@404	336
kpeter@408	337	/// \brief Sets the target node.
alpar@404	338	///
kpeter@408	339	/// Sets the target node.
kpeter@408	340	/// \return <tt>(*this)</tt>
kpeter@408	341	Preflow& target(const Node& node) {
kpeter@408	342	_target = node;
kpeter@408	343	return *this;
kpeter@408	344	}
kpeter@408	345
kpeter@408	346	/// \brief Sets the elevator used by algorithm.
kpeter@408	347	///
kpeter@408	348	/// Sets the elevator used by algorithm.
kpeter@408	349	/// If you don't use this function before calling \ref run() or
kpeter@408	350	/// \ref init(), an instance will be allocated automatically.
kpeter@408	351	/// The destructor deallocates this automatically allocated elevator,
kpeter@408	352	/// of course.
kpeter@408	353	/// \return <tt>(*this)</tt>
alpar@404	354	Preflow& elevator(Elevator& elevator) {
alpar@404	355	if (_local_level) {
alpar@404	356	delete _level;
alpar@404	357	_local_level = false;
alpar@404	358	}
alpar@404	359	_level = &elevator;
alpar@404	360	return *this;
alpar@404	361	}
alpar@404	362
kpeter@408	363	/// \brief Returns a const reference to the elevator.
alpar@404	364	///
kpeter@408	365	/// Returns a const reference to the elevator.
kpeter@408	366	///
kpeter@408	367	/// \pre Either \ref run() or \ref init() must be called before
kpeter@408	368	/// using this function.
kpeter@437	369	const Elevator& elevator() const {
alpar@404	370	return *_level;
alpar@404	371	}
alpar@404	372
alpar@404	373	/// \brief Sets the tolerance used by algorithm.
alpar@404	374	///
alpar@404	375	/// Sets the tolerance used by algorithm.
alpar@404	376	Preflow& tolerance(const Tolerance& tolerance) const {
alpar@404	377	_tolerance = tolerance;
alpar@404	378	return *this;
alpar@404	379	}
alpar@404	380
kpeter@408	381	/// \brief Returns a const reference to the tolerance.
alpar@404	382	///
kpeter@408	383	/// Returns a const reference to the tolerance.
alpar@404	384	const Tolerance& tolerance() const {
alpar@404	385	return tolerance;
alpar@404	386	}
alpar@404	387
kpeter@408	388	/// \name Execution Control
kpeter@408	389	/// The simplest way to execute the preflow algorithm is to use
kpeter@408	390	/// \ref run() or \ref runMinCut().\n
kpeter@408	391	/// If you need more control on the initial solution or the execution,
kpeter@408	392	/// first you have to call one of the \ref init() functions, then
kpeter@408	393	/// \ref startFirstPhase() and if you need it \ref startSecondPhase().
alpar@404	394
alpar@404	395	///@{
alpar@404	396
alpar@404	397	/// \brief Initializes the internal data structures.
alpar@404	398	///
kpeter@408	399	/// Initializes the internal data structures and sets the initial
kpeter@408	400	/// flow to zero on each arc.
alpar@404	401	void init() {
alpar@404	402	createStructures();
alpar@404	403
alpar@404	404	_phase = true;
alpar@404	405	for (NodeIt n(_graph); n != INVALID; ++n) {
alpar@404	406	_excess->set(n, 0);
alpar@404	407	}
alpar@404	408
alpar@404	409	for (ArcIt e(_graph); e != INVALID; ++e) {
alpar@404	410	_flow->set(e, 0);
alpar@404	411	}
alpar@404	412
alpar@404	413	typename Digraph::template NodeMap<bool> reached(_graph, false);
alpar@404	414
alpar@404	415	_level->initStart();
alpar@404	416	_level->initAddItem(_target);
alpar@404	417
alpar@404	418	std::vector<Node> queue;
alpar@404	419	reached.set(_source, true);
alpar@404	420
alpar@404	421	queue.push_back(_target);
alpar@404	422	reached.set(_target, true);
alpar@404	423	while (!queue.empty()) {
alpar@404	424	_level->initNewLevel();
alpar@404	425	std::vector<Node> nqueue;
alpar@404	426	for (int i = 0; i < int(queue.size()); ++i) {
alpar@404	427	Node n = queue[i];
alpar@404	428	for (InArcIt e(_graph, n); e != INVALID; ++e) {
alpar@404	429	Node u = _graph.source(e);
alpar@404	430	if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
alpar@404	431	reached.set(u, true);
alpar@404	432	_level->initAddItem(u);
alpar@404	433	nqueue.push_back(u);
alpar@404	434	}
alpar@404	435	}
alpar@404	436	}
alpar@404	437	queue.swap(nqueue);
alpar@404	438	}
alpar@404	439	_level->initFinish();
alpar@404	440
alpar@404	441	for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
alpar@404	442	if (_tolerance.positive((*_capacity)[e])) {
alpar@404	443	Node u = _graph.target(e);
alpar@404	444	if ((*_level)[u] == _level->maxLevel()) continue;
alpar@404	445	_flow->set(e, (*_capacity)[e]);
alpar@404	446	_excess->set(u, (_excess)[u] + (_capacity)[e]);
alpar@404	447	if (u != _target && !_level->active(u)) {
alpar@404	448	_level->activate(u);
alpar@404	449	}
alpar@404	450	}
alpar@404	451	}
alpar@404	452	}
alpar@404	453
kpeter@408	454	/// \brief Initializes the internal data structures using the
kpeter@408	455	/// given flow map.
alpar@404	456	///
alpar@404	457	/// Initializes the internal data structures and sets the initial
alpar@404	458	/// flow to the given \c flowMap. The \c flowMap should contain a
kpeter@408	459	/// flow or at least a preflow, i.e. at each node excluding the
kpeter@408	460	/// source node the incoming flow should greater or equal to the
alpar@404	461	/// outgoing flow.
kpeter@408	462	/// \return \c false if the given \c flowMap is not a preflow.
alpar@404	463	template <typename FlowMap>
kpeter@407	464	bool init(const FlowMap& flowMap) {
alpar@404	465	createStructures();
alpar@404	466
alpar@404	467	for (ArcIt e(_graph); e != INVALID; ++e) {
alpar@404	468	_flow->set(e, flowMap[e]);
alpar@404	469	}
alpar@404	470
alpar@404	471	for (NodeIt n(_graph); n != INVALID; ++n) {
alpar@404	472	Value excess = 0;
alpar@404	473	for (InArcIt e(_graph, n); e != INVALID; ++e) {
alpar@404	474	excess += (*_flow)[e];
alpar@404	475	}
alpar@404	476	for (OutArcIt e(_graph, n); e != INVALID; ++e) {
alpar@404	477	excess -= (*_flow)[e];
alpar@404	478	}
alpar@404	479	if (excess < 0 && n != _source) return false;
alpar@404	480	_excess->set(n, excess);
alpar@404	481	}
alpar@404	482
alpar@404	483	typename Digraph::template NodeMap<bool> reached(_graph, false);
alpar@404	484
alpar@404	485	_level->initStart();
alpar@404	486	_level->initAddItem(_target);
alpar@404	487
alpar@404	488	std::vector<Node> queue;
alpar@404	489	reached.set(_source, true);
alpar@404	490
alpar@404	491	queue.push_back(_target);
alpar@404	492	reached.set(_target, true);
alpar@404	493	while (!queue.empty()) {
alpar@404	494	_level->initNewLevel();
alpar@404	495	std::vector<Node> nqueue;
alpar@404	496	for (int i = 0; i < int(queue.size()); ++i) {
alpar@404	497	Node n = queue[i];
alpar@404	498	for (InArcIt e(_graph, n); e != INVALID; ++e) {
alpar@404	499	Node u = _graph.source(e);
alpar@404	500	if (!reached[u] &&
alpar@404	501	_tolerance.positive((_capacity)[e] - (_flow)[e])) {
alpar@404	502	reached.set(u, true);
alpar@404	503	_level->initAddItem(u);
alpar@404	504	nqueue.push_back(u);
alpar@404	505	}
alpar@404	506	}
alpar@404	507	for (OutArcIt e(_graph, n); e != INVALID; ++e) {
alpar@404	508	Node v = _graph.target(e);
alpar@404	509	if (!reached[v] && _tolerance.positive((*_flow)[e])) {
alpar@404	510	reached.set(v, true);
alpar@404	511	_level->initAddItem(v);
alpar@404	512	nqueue.push_back(v);
alpar@404	513	}
alpar@404	514	}
alpar@404	515	}
alpar@404	516	queue.swap(nqueue);
alpar@404	517	}
alpar@404	518	_level->initFinish();
alpar@404	519
alpar@404	520	for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
alpar@404	521	Value rem = (_capacity)[e] - (_flow)[e];
alpar@404	522	if (_tolerance.positive(rem)) {
alpar@404	523	Node u = _graph.target(e);
alpar@404	524	if ((*_level)[u] == _level->maxLevel()) continue;
alpar@404	525	_flow->set(e, (*_capacity)[e]);
alpar@404	526	_excess->set(u, (*_excess)[u] + rem);
alpar@404	527	if (u != _target && !_level->active(u)) {
alpar@404	528	_level->activate(u);
alpar@404	529	}
alpar@404	530	}
alpar@404	531	}
alpar@404	532	for (InArcIt e(_graph, _source); e != INVALID; ++e) {
alpar@404	533	Value rem = (*_flow)[e];
alpar@404	534	if (_tolerance.positive(rem)) {
alpar@404	535	Node v = _graph.source(e);
alpar@404	536	if ((*_level)[v] == _level->maxLevel()) continue;
alpar@404	537	_flow->set(e, 0);
alpar@404	538	_excess->set(v, (*_excess)[v] + rem);
alpar@404	539	if (v != _target && !_level->active(v)) {
alpar@404	540	_level->activate(v);
alpar@404	541	}
alpar@404	542	}
alpar@404	543	}
alpar@404	544	return true;
alpar@404	545	}
alpar@404	546
alpar@404	547	/// \brief Starts the first phase of the preflow algorithm.
alpar@404	548	///
alpar@404	549	/// The preflow algorithm consists of two phases, this method runs
alpar@404	550	/// the first phase. After the first phase the maximum flow value
alpar@404	551	/// and a minimum value cut can already be computed, although a
alpar@404	552	/// maximum flow is not yet obtained. So after calling this method
alpar@404	553	/// \ref flowValue() returns the value of a maximum flow and \ref
alpar@404	554	/// minCut() returns a minimum cut.
kpeter@408	555	/// \pre One of the \ref init() functions must be called before
kpeter@408	556	/// using this function.
alpar@404	557	void startFirstPhase() {
alpar@404	558	_phase = true;
alpar@404	559
alpar@404	560	Node n = _level->highestActive();
alpar@404	561	int level = _level->highestActiveLevel();
alpar@404	562	while (n != INVALID) {
alpar@404	563	int num = _node_num;
alpar@404	564
alpar@404	565	while (num > 0 && n != INVALID) {
alpar@404	566	Value excess = (*_excess)[n];
alpar@404	567	int new_level = _level->maxLevel();
alpar@404	568
alpar@404	569	for (OutArcIt e(_graph, n); e != INVALID; ++e) {
alpar@404	570	Value rem = (_capacity)[e] - (_flow)[e];
alpar@404	571	if (!_tolerance.positive(rem)) continue;
alpar@404	572	Node v = _graph.target(e);
alpar@404	573	if ((*_level)[v] < level) {
alpar@404	574	if (!_level->active(v) && v != _target) {
alpar@404	575	_level->activate(v);
alpar@404	576	}
alpar@404	577	if (!_tolerance.less(rem, excess)) {
alpar@404	578	_flow->set(e, (*_flow)[e] + excess);
alpar@404	579	_excess->set(v, (*_excess)[v] + excess);
alpar@404	580	excess = 0;
alpar@404	581	goto no_more_push_1;
alpar@404	582	} else {
alpar@404	583	excess -= rem;
alpar@404	584	_excess->set(v, (*_excess)[v] + rem);
alpar@404	585	_flow->set(e, (*_capacity)[e]);
alpar@404	586	}
alpar@404	587	} else if (new_level > (*_level)[v]) {
alpar@404	588	new_level = (*_level)[v];
alpar@404	589	}
alpar@404	590	}
alpar@404	591
alpar@404	592	for (InArcIt e(_graph, n); e != INVALID; ++e) {
alpar@404	593	Value rem = (*_flow)[e];
alpar@404	594	if (!_tolerance.positive(rem)) continue;
alpar@404	595	Node v = _graph.source(e);
alpar@404	596	if ((*_level)[v] < level) {
alpar@404	597	if (!_level->active(v) && v != _target) {
alpar@404	598	_level->activate(v);
alpar@404	599	}
alpar@404	600	if (!_tolerance.less(rem, excess)) {
alpar@404	601	_flow->set(e, (*_flow)[e] - excess);
alpar@404	602	_excess->set(v, (*_excess)[v] + excess);
alpar@404	603	excess = 0;
alpar@404	604	goto no_more_push_1;
alpar@404	605	} else {
alpar@404	606	excess -= rem;
alpar@404	607	_excess->set(v, (*_excess)[v] + rem);
alpar@404	608	_flow->set(e, 0);
alpar@404	609	}
alpar@404	610	} else if (new_level > (*_level)[v]) {
alpar@404	611	new_level = (*_level)[v];
alpar@404	612	}
alpar@404	613	}
alpar@404	614
alpar@404	615	no_more_push_1:
alpar@404	616
alpar@404	617	_excess->set(n, excess);
alpar@404	618
alpar@404	619	if (excess != 0) {
alpar@404	620	if (new_level + 1 < _level->maxLevel()) {
alpar@404	621	_level->liftHighestActive(new_level + 1);
alpar@404	622	} else {
alpar@404	623	_level->liftHighestActiveToTop();
alpar@404	624	}
alpar@404	625	if (_level->emptyLevel(level)) {
alpar@404	626	_level->liftToTop(level);
alpar@404	627	}
alpar@404	628	} else {
alpar@404	629	_level->deactivate(n);
alpar@404	630	}
alpar@404	631
alpar@404	632	n = _level->highestActive();
alpar@404	633	level = _level->highestActiveLevel();
alpar@404	634	--num;
alpar@404	635	}
alpar@404	636
alpar@404	637	num = _node_num * 20;
alpar@404	638	while (num > 0 && n != INVALID) {
alpar@404	639	Value excess = (*_excess)[n];
alpar@404	640	int new_level = _level->maxLevel();
alpar@404	641
alpar@404	642	for (OutArcIt e(_graph, n); e != INVALID; ++e) {
alpar@404	643	Value rem = (_capacity)[e] - (_flow)[e];
alpar@404	644	if (!_tolerance.positive(rem)) continue;
alpar@404	645	Node v = _graph.target(e);
alpar@404	646	if ((*_level)[v] < level) {
alpar@404	647	if (!_level->active(v) && v != _target) {
alpar@404	648	_level->activate(v);
alpar@404	649	}
alpar@404	650	if (!_tolerance.less(rem, excess)) {
alpar@404	651	_flow->set(e, (*_flow)[e] + excess);
alpar@404	652	_excess->set(v, (*_excess)[v] + excess);
alpar@404	653	excess = 0;
alpar@404	654	goto no_more_push_2;
alpar@404	655	} else {
alpar@404	656	excess -= rem;
alpar@404	657	_excess->set(v, (*_excess)[v] + rem);
alpar@404	658	_flow->set(e, (*_capacity)[e]);
alpar@404	659	}
alpar@404	660	} else if (new_level > (*_level)[v]) {
alpar@404	661	new_level = (*_level)[v];
alpar@404	662	}
alpar@404	663	}
alpar@404	664
alpar@404	665	for (InArcIt e(_graph, n); e != INVALID; ++e) {
alpar@404	666	Value rem = (*_flow)[e];
alpar@404	667	if (!_tolerance.positive(rem)) continue;
alpar@404	668	Node v = _graph.source(e);
alpar@404	669	if ((*_level)[v] < level) {
alpar@404	670	if (!_level->active(v) && v != _target) {
alpar@404	671	_level->activate(v);
alpar@404	672	}
alpar@404	673	if (!_tolerance.less(rem, excess)) {
alpar@404	674	_flow->set(e, (*_flow)[e] - excess);
alpar@404	675	_excess->set(v, (*_excess)[v] + excess);
alpar@404	676	excess = 0;
alpar@404	677	goto no_more_push_2;
alpar@404	678	} else {
alpar@404	679	excess -= rem;
alpar@404	680	_excess->set(v, (*_excess)[v] + rem);
alpar@404	681	_flow->set(e, 0);
alpar@404	682	}
alpar@404	683	} else if (new_level > (*_level)[v]) {
alpar@404	684	new_level = (*_level)[v];
alpar@404	685	}
alpar@404	686	}
alpar@404	687
alpar@404	688	no_more_push_2:
alpar@404	689
alpar@404	690	_excess->set(n, excess);
alpar@404	691
alpar@404	692	if (excess != 0) {
alpar@404	693	if (new_level + 1 < _level->maxLevel()) {
alpar@404	694	_level->liftActiveOn(level, new_level + 1);
alpar@404	695	} else {
alpar@404	696	_level->liftActiveToTop(level);
alpar@404	697	}
alpar@404	698	if (_level->emptyLevel(level)) {
alpar@404	699	_level->liftToTop(level);
alpar@404	700	}
alpar@404	701	} else {
alpar@404	702	_level->deactivate(n);
alpar@404	703	}
alpar@404	704
alpar@404	705	while (level >= 0 && _level->activeFree(level)) {
alpar@404	706	--level;
alpar@404	707	}
alpar@404	708	if (level == -1) {
alpar@404	709	n = _level->highestActive();
alpar@404	710	level = _level->highestActiveLevel();
alpar@404	711	} else {
alpar@404	712	n = _level->activeOn(level);
alpar@404	713	}
alpar@404	714	--num;
alpar@404	715	}
alpar@404	716	}
alpar@404	717	}
alpar@404	718
alpar@404	719	/// \brief Starts the second phase of the preflow algorithm.
alpar@404	720	///
alpar@404	721	/// The preflow algorithm consists of two phases, this method runs
kpeter@408	722	/// the second phase. After calling one of the \ref init() functions
kpeter@408	723	/// and \ref startFirstPhase() and then \ref startSecondPhase(),
kpeter@408	724	/// \ref flowMap() returns a maximum flow, \ref flowValue() returns the
alpar@404	725	/// value of a maximum flow, \ref minCut() returns a minimum cut
kpeter@408	726	/// \pre One of the \ref init() functions and \ref startFirstPhase()
kpeter@408	727	/// must be called before using this function.
alpar@404	728	void startSecondPhase() {
alpar@404	729	_phase = false;
alpar@404	730
alpar@404	731	typename Digraph::template NodeMap<bool> reached(_graph);
alpar@404	732	for (NodeIt n(_graph); n != INVALID; ++n) {
alpar@404	733	reached.set(n, (*_level)[n] < _level->maxLevel());
alpar@404	734	}
alpar@404	735
alpar@404	736	_level->initStart();
alpar@404	737	_level->initAddItem(_source);
alpar@404	738
alpar@404	739	std::vector<Node> queue;
alpar@404	740	queue.push_back(_source);
alpar@404	741	reached.set(_source, true);
alpar@404	742
alpar@404	743	while (!queue.empty()) {
alpar@404	744	_level->initNewLevel();
alpar@404	745	std::vector<Node> nqueue;
alpar@404	746	for (int i = 0; i < int(queue.size()); ++i) {
alpar@404	747	Node n = queue[i];
alpar@404	748	for (OutArcIt e(_graph, n); e != INVALID; ++e) {
alpar@404	749	Node v = _graph.target(e);
alpar@404	750	if (!reached[v] && _tolerance.positive((*_flow)[e])) {
alpar@404	751	reached.set(v, true);
alpar@404	752	_level->initAddItem(v);
alpar@404	753	nqueue.push_back(v);
alpar@404	754	}
alpar@404	755	}
alpar@404	756	for (InArcIt e(_graph, n); e != INVALID; ++e) {
alpar@404	757	Node u = _graph.source(e);
alpar@404	758	if (!reached[u] &&
alpar@404	759	_tolerance.positive((_capacity)[e] - (_flow)[e])) {
alpar@404	760	reached.set(u, true);
alpar@404	761	_level->initAddItem(u);
alpar@404	762	nqueue.push_back(u);
alpar@404	763	}
alpar@404	764	}
alpar@404	765	}
alpar@404	766	queue.swap(nqueue);
alpar@404	767	}
alpar@404	768	_level->initFinish();
alpar@404	769
alpar@404	770	for (NodeIt n(_graph); n != INVALID; ++n) {
alpar@404	771	if (!reached[n]) {
alpar@404	772	_level->dirtyTopButOne(n);
alpar@404	773	} else if ((*_excess)[n] > 0 && _target != n) {
alpar@404	774	_level->activate(n);
alpar@404	775	}
alpar@404	776	}
alpar@404	777
alpar@404	778	Node n;
alpar@404	779	while ((n = _level->highestActive()) != INVALID) {
alpar@404	780	Value excess = (*_excess)[n];
alpar@404	781	int level = _level->highestActiveLevel();
alpar@404	782	int new_level = _level->maxLevel();
alpar@404	783
alpar@404	784	for (OutArcIt e(_graph, n); e != INVALID; ++e) {
alpar@404	785	Value rem = (_capacity)[e] - (_flow)[e];
alpar@404	786	if (!_tolerance.positive(rem)) continue;
alpar@404	787	Node v = _graph.target(e);
alpar@404	788	if ((*_level)[v] < level) {
alpar@404	789	if (!_level->active(v) && v != _source) {
alpar@404	790	_level->activate(v);
alpar@404	791	}
alpar@404	792	if (!_tolerance.less(rem, excess)) {
alpar@404	793	_flow->set(e, (*_flow)[e] + excess);
alpar@404	794	_excess->set(v, (*_excess)[v] + excess);
alpar@404	795	excess = 0;
alpar@404	796	goto no_more_push;
alpar@404	797	} else {
alpar@404	798	excess -= rem;
alpar@404	799	_excess->set(v, (*_excess)[v] + rem);
alpar@404	800	_flow->set(e, (*_capacity)[e]);
alpar@404	801	}
alpar@404	802	} else if (new_level > (*_level)[v]) {
alpar@404	803	new_level = (*_level)[v];
alpar@404	804	}
alpar@404	805	}
alpar@404	806
alpar@404	807	for (InArcIt e(_graph, n); e != INVALID; ++e) {
alpar@404	808	Value rem = (*_flow)[e];
alpar@404	809	if (!_tolerance.positive(rem)) continue;
alpar@404	810	Node v = _graph.source(e);
alpar@404	811	if ((*_level)[v] < level) {
alpar@404	812	if (!_level->active(v) && v != _source) {
alpar@404	813	_level->activate(v);
alpar@404	814	}
alpar@404	815	if (!_tolerance.less(rem, excess)) {
alpar@404	816	_flow->set(e, (*_flow)[e] - excess);
alpar@404	817	_excess->set(v, (*_excess)[v] + excess);
alpar@404	818	excess = 0;
alpar@404	819	goto no_more_push;
alpar@404	820	} else {
alpar@404	821	excess -= rem;
alpar@404	822	_excess->set(v, (*_excess)[v] + rem);
alpar@404	823	_flow->set(e, 0);
alpar@404	824	}
alpar@404	825	} else if (new_level > (*_level)[v]) {
alpar@404	826	new_level = (*_level)[v];
alpar@404	827	}
alpar@404	828	}
alpar@404	829
alpar@404	830	no_more_push:
alpar@404	831
alpar@404	832	_excess->set(n, excess);
alpar@404	833
alpar@404	834	if (excess != 0) {
alpar@404	835	if (new_level + 1 < _level->maxLevel()) {
alpar@404	836	_level->liftHighestActive(new_level + 1);
alpar@404	837	} else {
alpar@404	838	// Calculation error
alpar@404	839	_level->liftHighestActiveToTop();
alpar@404	840	}
alpar@404	841	if (_level->emptyLevel(level)) {
alpar@404	842	// Calculation error
alpar@404	843	_level->liftToTop(level);
alpar@404	844	}
alpar@404	845	} else {
alpar@404	846	_level->deactivate(n);
alpar@404	847	}
alpar@404	848
alpar@404	849	}
alpar@404	850	}
alpar@404	851
alpar@404	852	/// \brief Runs the preflow algorithm.
alpar@404	853	///
alpar@404	854	/// Runs the preflow algorithm.
alpar@404	855	/// \note pf.run() is just a shortcut of the following code.
alpar@404	856	/// \code
alpar@404	857	/// pf.init();
alpar@404	858	/// pf.startFirstPhase();
alpar@404	859	/// pf.startSecondPhase();
alpar@404	860	/// \endcode
alpar@404	861	void run() {
alpar@404	862	init();
alpar@404	863	startFirstPhase();
alpar@404	864	startSecondPhase();
alpar@404	865	}
alpar@404	866
alpar@404	867	/// \brief Runs the preflow algorithm to compute the minimum cut.
alpar@404	868	///
alpar@404	869	/// Runs the preflow algorithm to compute the minimum cut.
alpar@404	870	/// \note pf.runMinCut() is just a shortcut of the following code.
alpar@404	871	/// \code
alpar@404	872	/// pf.init();
alpar@404	873	/// pf.startFirstPhase();
alpar@404	874	/// \endcode
alpar@404	875	void runMinCut() {
alpar@404	876	init();
alpar@404	877	startFirstPhase();
alpar@404	878	}
alpar@404	879
alpar@404	880	/// @}
alpar@404	881
alpar@404	882	/// \name Query Functions
kpeter@408	883	/// The results of the preflow algorithm can be obtained using these
alpar@404	884	/// functions.\n
kpeter@408	885	/// Either one of the \ref run() "run*()" functions or one of the
kpeter@408	886	/// \ref startFirstPhase() "start*()" functions should be called
kpeter@408	887	/// before using them.
alpar@404	888
alpar@404	889	///@{
alpar@404	890
alpar@404	891	/// \brief Returns the value of the maximum flow.
alpar@404	892	///
alpar@404	893	/// Returns the value of the maximum flow by returning the excess
kpeter@408	894	/// of the target node. This value equals to the value of
kpeter@408	895	/// the maximum flow already after the first phase of the algorithm.
kpeter@408	896	///
kpeter@408	897	/// \pre Either \ref run() or \ref init() must be called before
kpeter@408	898	/// using this function.
alpar@404	899	Value flowValue() const {
alpar@404	900	return (*_excess)[_target];
alpar@404	901	}
alpar@404	902
kpeter@408	903	/// \brief Returns the flow on the given arc.
alpar@404	904	///
kpeter@408	905	/// Returns the flow on the given arc. This method can
kpeter@408	906	/// be called after the second phase of the algorithm.
kpeter@408	907	///
kpeter@408	908	/// \pre Either \ref run() or \ref init() must be called before
kpeter@408	909	/// using this function.
kpeter@408	910	Value flow(const Arc& arc) const {
kpeter@408	911	return (*_flow)[arc];
kpeter@408	912	}
kpeter@408	913
kpeter@408	914	/// \brief Returns a const reference to the flow map.
kpeter@408	915	///
kpeter@408	916	/// Returns a const reference to the arc map storing the found flow.
kpeter@408	917	/// This method can be called after the second phase of the algorithm.
kpeter@408	918	///
kpeter@408	919	/// \pre Either \ref run() or \ref init() must be called before
kpeter@408	920	/// using this function.
kpeter@437	921	const FlowMap& flowMap() const {
kpeter@408	922	return *_flow;
kpeter@408	923	}
kpeter@408	924
kpeter@408	925	/// \brief Returns \c true when the node is on the source side of the
kpeter@408	926	/// minimum cut.
kpeter@408	927	///
kpeter@408	928	/// Returns true when the node is on the source side of the found
kpeter@408	929	/// minimum cut. This method can be called both after running \ref
alpar@404	930	/// startFirstPhase() and \ref startSecondPhase().
kpeter@408	931	///
kpeter@408	932	/// \pre Either \ref run() or \ref init() must be called before
kpeter@408	933	/// using this function.
alpar@404	934	bool minCut(const Node& node) const {
alpar@404	935	return ((*_level)[node] == _level->maxLevel()) == _phase;
alpar@404	936	}
alpar@404	937
kpeter@408	938	/// \brief Gives back a minimum value cut.
alpar@404	939	///
kpeter@408	940	/// Sets \c cutMap to the characteristic vector of a minimum value
kpeter@408	941	/// cut. \c cutMap should be a \ref concepts::WriteMap "writable"
kpeter@408	942	/// node map with \c bool (or convertible) value type.
kpeter@408	943	///
kpeter@408	944	/// This method can be called both after running \ref startFirstPhase()
kpeter@408	945	/// and \ref startSecondPhase(). The result after the second phase
kpeter@408	946	/// could be slightly different if inexact computation is used.
kpeter@408	947	///
kpeter@408	948	/// \note This function calls \ref minCut() for each node, so it runs in
kpeter@408	949	/// \f$O(n)\f$ time.
kpeter@408	950	///
kpeter@408	951	/// \pre Either \ref run() or \ref init() must be called before
kpeter@408	952	/// using this function.
alpar@404	953	template <typename CutMap>
alpar@404	954	void minCutMap(CutMap& cutMap) const {
alpar@404	955	for (NodeIt n(_graph); n != INVALID; ++n) {
alpar@404	956	cutMap.set(n, minCut(n));
alpar@404	957	}
alpar@404	958	}
alpar@404	959
alpar@404	960	/// @}
alpar@404	961	};
alpar@404	962	}
alpar@404	963
alpar@404	964	#endif

author	Peter Kovacs <kpeter@inf.elte.hu>
	Mon, 12 Jan 2009 08:18:04 +0100
changeset 488	2b5496c62ccd
parent 437	6a2a33ad261b
child 525	9605e051942f
permissions	-rw-r--r--