lemon-0.x: src/work/jacint/preflow_max_flow.h@cf7b01232d86 (annotated)

jacint@109	1	// -- C++ --
jacint@109	2	/*
jacint@109	3	preflow_max_flow.h
jacint@109	4	by jacint.
jacint@109	5	Runs the first phase of preflow.h
jacint@109	6
jacint@109	7	The constructor runs the algorithm.
jacint@109	8
jacint@109	9	Members:
jacint@109	10
jacint@109	11	T maxFlow() : returns the value of a maximum flow
jacint@109	12
jacint@109	13	CutMap minCut() : returns the characteristic vector of a min cut.
jacint@109	14	*/
jacint@109	15
jacint@109	16	#ifndef PREFLOW_MAX_FLOW_H
jacint@109	17	#define PREFLOW_MAX_FLOW_H
jacint@109	18
jacint@109	19	#define H0 20
jacint@109	20	#define H1 1
jacint@109	21
jacint@109	22	#include <vector>
jacint@109	23	#include <queue>
jacint@109	24
jacint@109	25	namespace hugo {
jacint@109	26
jacint@109	27	template <typename Graph, typename T,
jacint@109	28	typename FlowMap=typename Graph::EdgeMap<T>,
jacint@109	29	typename CapMap=typename Graph::EdgeMap<T>,
jacint@109	30	typename CutMap=typename Graph::NodeMap<bool> >
jacint@109	31	class preflow_max_flow {
jacint@109	32
jacint@109	33	typedef typename Graph::NodeIt NodeIt;
jacint@109	34	typedef typename Graph::EdgeIt EdgeIt;
jacint@109	35	typedef typename Graph::EachNodeIt EachNodeIt;
jacint@109	36	typedef typename Graph::OutEdgeIt OutEdgeIt;
jacint@109	37	typedef typename Graph::InEdgeIt InEdgeIt;
jacint@109	38
jacint@109	39	Graph& G;
jacint@109	40	NodeIt s;
jacint@109	41	NodeIt t;
jacint@109	42	FlowMap flow;
jacint@109	43	CapMap& capacity;
jacint@109	44	CutMap cut;
jacint@109	45	T value;
jacint@109	46
jacint@109	47	public:
jacint@109	48
jacint@109	49	preflow_max_flow(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity ) :
jacint@109	50	G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity), cut(_G, false)
jacint@109	51	{
jacint@109	52
jacint@109	53	int n=G.nodeNum();
jacint@109	54	int heur0=(int)(H0*n); //time while running 'bound decrease'
jacint@109	55	int heur1=(int)(H1*n); //time while running 'highest label'
jacint@109	56	int heur=heur1; //starting time interval (#of relabels)
jacint@109	57	bool what_heur=1;
jacint@109	58	/*
jacint@109	59	what_heur is 0 in case 'bound decrease'
jacint@109	60	and 1 in case 'highest label'
jacint@109	61	*/
jacint@109	62	bool end=false;
jacint@109	63	/*
jacint@109	64	Needed for 'bound decrease', 'true'
jacint@109	65	means no active nodes are above bound b.
jacint@109	66	*/
jacint@109	67	int relabel=0;
jacint@109	68	int k=n-2; //bound on the highest level under n containing a node
jacint@109	69	int b=k; //bound on the highest level under n of an active node
jacint@109	70
jacint@109	71	typename Graph::NodeMap<int> level(G,n);
jacint@109	72	typename Graph::NodeMap<T> excess(G);
jacint@109	73
jacint@109	74	std::vector<NodeIt> active(n);
jacint@109	75	typename Graph::NodeMap<NodeIt> next(G);
jacint@109	76	//Stack of the active nodes in level i < n.
jacint@109	77	//We use it in both phases.
jacint@109	78
jacint@109	79	typename Graph::NodeMap<NodeIt> left(G);
jacint@109	80	typename Graph::NodeMap<NodeIt> right(G);
jacint@109	81	std::vector<NodeIt> level_list(n);
jacint@109	82	/*
jacint@109	83	List of the nodes in level i<n.
jacint@109	84	*/
jacint@109	85
jacint@109	86	/Reverse_bfs from t, to find the starting level./
jacint@109	87	level.set(t,0);
jacint@109	88	std::queue<NodeIt> bfs_queue;
jacint@109	89	bfs_queue.push(t);
jacint@109	90
jacint@109	91	while (!bfs_queue.empty()) {
jacint@109	92
jacint@109	93	NodeIt v=bfs_queue.front();
jacint@109	94	bfs_queue.pop();
jacint@109	95	int l=level.get(v)+1;
jacint@109	96
jacint@109	97	for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
jacint@109	98	NodeIt w=G.tail(e);
jacint@109	99	if ( level.get(w) == n && w != s ) {
jacint@109	100	bfs_queue.push(w);
jacint@109	101	NodeIt first=level_list[l];
jacint@109	102	if ( first != 0 ) left.set(first,w);
jacint@109	103	right.set(w,first);
jacint@109	104	level_list[l]=w;
jacint@109	105	level.set(w, l);
jacint@109	106	}
jacint@109	107	}
jacint@109	108	}
jacint@109	109
jacint@109	110	level.set(s,n);
jacint@109	111
jacint@109	112
jacint@109	113	/* Starting flow. It is everywhere 0 at the moment. */
jacint@109	114	for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e)
jacint@109	115	{
jacint@109	116	T c=capacity.get(e);
jacint@109	117	if ( c == 0 ) continue;
jacint@109	118	NodeIt w=G.head(e);
jacint@109	119	if ( level.get(w) < n ) {
jacint@109	120	if ( excess.get(w) == 0 && w!=t ) {
jacint@109	121	next.set(w,active[level.get(w)]);
jacint@109	122	active[level.get(w)]=w;
jacint@109	123	}
jacint@109	124	flow.set(e, c);
jacint@109	125	excess.set(w, excess.get(w)+c);
jacint@109	126	}
jacint@109	127	}
jacint@109	128
jacint@109	129	/*
jacint@109	130	End of preprocessing
jacint@109	131	*/
jacint@109	132
jacint@109	133
jacint@109	134
jacint@109	135	/*
jacint@109	136	Push/relabel on the highest level active nodes.
jacint@109	137	*/
jacint@109	138	while ( true ) {
jacint@109	139
jacint@109	140	if ( b == 0 ) {
jacint@109	141	if ( !what_heur && !end && k > 0 ) {
jacint@109	142	b=k;
jacint@109	143	end=true;
jacint@109	144	} else break;
jacint@109	145	}
jacint@109	146
jacint@109	147
jacint@109	148	if ( active[b] == 0 ) --b;
jacint@109	149	else {
jacint@109	150	end=false;
jacint@109	151
jacint@109	152	NodeIt w=active[b];
jacint@109	153	active[b]=next.get(w);
jacint@109	154	int lev=level.get(w);
jacint@109	155	T exc=excess.get(w);
jacint@109	156	int newlevel=n; //bound on the next level of w
jacint@109	157
jacint@109	158	for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
jacint@109	159
jacint@109	160	if ( flow.get(e) == capacity.get(e) ) continue;
jacint@109	161	NodeIt v=G.head(e);
jacint@109	162	//e=wv
jacint@109	163
jacint@109	164	if( lev > level.get(v) ) {
jacint@109	165	/Push is allowed now/
jacint@109	166
jacint@109	167	if ( excess.get(v)==0 && v!=t && v!=s ) {
jacint@109	168	int lev_v=level.get(v);
jacint@109	169	next.set(v,active[lev_v]);
jacint@109	170	active[lev_v]=v;
jacint@109	171	}
jacint@109	172
jacint@109	173	T cap=capacity.get(e);
jacint@109	174	T flo=flow.get(e);
jacint@109	175	T remcap=cap-flo;
jacint@109	176
jacint@109	177	if ( remcap >= exc ) {
jacint@109	178	/A nonsaturating push./
jacint@109	179
jacint@109	180	flow.set(e, flo+exc);
jacint@109	181	excess.set(v, excess.get(v)+exc);
jacint@109	182	exc=0;
jacint@109	183	break;
jacint@109	184
jacint@109	185	} else {
jacint@109	186	/A saturating push./
jacint@109	187
jacint@109	188	flow.set(e, cap);
jacint@109	189	excess.set(v, excess.get(v)+remcap);
jacint@109	190	exc-=remcap;
jacint@109	191	}
jacint@109	192	} else if ( newlevel > level.get(v) ){
jacint@109	193	newlevel = level.get(v);
jacint@109	194	}
jacint@109	195
jacint@109	196	} //for out edges wv
jacint@109	197
jacint@109	198
jacint@109	199	if ( exc > 0 ) {
jacint@109	200	for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
jacint@109	201
jacint@109	202	if( flow.get(e) == 0 ) continue;
jacint@109	203	NodeIt v=G.tail(e);
jacint@109	204	//e=vw
jacint@109	205
jacint@109	206	if( lev > level.get(v) ) {
jacint@109	207	/Push is allowed now/
jacint@109	208
jacint@109	209	if ( excess.get(v)==0 && v!=t && v!=s ) {
jacint@109	210	int lev_v=level.get(v);
jacint@109	211	next.set(v,active[lev_v]);
jacint@109	212	active[lev_v]=v;
jacint@109	213	}
jacint@109	214
jacint@109	215	T flo=flow.get(e);
jacint@109	216
jacint@109	217	if ( flo >= exc ) {
jacint@109	218	/A nonsaturating push./
jacint@109	219
jacint@109	220	flow.set(e, flo-exc);
jacint@109	221	excess.set(v, excess.get(v)+exc);
jacint@109	222	exc=0;
jacint@109	223	break;
jacint@109	224	} else {
jacint@109	225	/A saturating push./
jacint@109	226
jacint@109	227	excess.set(v, excess.get(v)+flo);
jacint@109	228	exc-=flo;
jacint@109	229	flow.set(e,0);
jacint@109	230	}
jacint@109	231	} else if ( newlevel > level.get(v) ) {
jacint@109	232	newlevel = level.get(v);
jacint@109	233	}
jacint@109	234	} //for in edges vw
jacint@109	235
jacint@109	236	} // if w still has excess after the out edge for cycle
jacint@109	237
jacint@109	238	excess.set(w, exc);
jacint@109	239
jacint@109	240	/*
jacint@109	241	Relabel
jacint@109	242	*/
jacint@109	243
jacint@109	244
jacint@109	245	if ( exc > 0 ) {
jacint@109	246	//now 'lev' is the old level of w
jacint@109	247
jacint@109	248	//unlacing starts
jacint@109	249	NodeIt right_n=right.get(w);
jacint@109	250	NodeIt left_n=left.get(w);
jacint@109	251
jacint@109	252	if ( right_n != 0 ) {
jacint@109	253	if ( left_n != 0 ) {
jacint@109	254	right.set(left_n, right_n);
jacint@109	255	left.set(right_n, left_n);
jacint@109	256	} else {
jacint@109	257	level_list[lev]=right_n;
jacint@109	258	left.set(right_n, 0);
jacint@109	259	}
jacint@109	260	} else {
jacint@109	261	if ( left_n != 0 ) {
jacint@109	262	right.set(left_n, 0);
jacint@109	263	} else {
jacint@109	264	level_list[lev]=0;
jacint@109	265
jacint@109	266	}
jacint@109	267	}
jacint@109	268	//unlacing ends
jacint@109	269
jacint@109	270	//gapping starts
jacint@109	271	if ( level_list[lev]==0 ) {
jacint@109	272
jacint@109	273	for (int i=lev; i!=k ; ) {
jacint@109	274	NodeIt v=level_list[++i];
jacint@109	275	while ( v != 0 ) {
jacint@109	276	level.set(v,n);
jacint@109	277	v=right.get(v);
jacint@109	278	}
jacint@109	279	level_list[i]=0;
jacint@109	280	if ( !what_heur ) active[i]=0;
jacint@109	281	}
jacint@109	282
jacint@109	283	level.set(w,n);
jacint@109	284	b=lev-1;
jacint@109	285	k=b;
jacint@109	286	//gapping ends
jacint@109	287	} else {
jacint@109	288
jacint@109	289	if ( newlevel == n ) level.set(w,n);
jacint@109	290	else {
jacint@109	291	level.set(w,++newlevel);
jacint@109	292	next.set(w,active[newlevel]);
jacint@109	293	active[newlevel]=w;
jacint@109	294	if ( what_heur ) b=newlevel;
jacint@109	295	if ( k < newlevel ) ++k;
jacint@109	296	NodeIt first=level_list[newlevel];
jacint@109	297	if ( first != 0 ) left.set(first,w);
jacint@109	298	right.set(w,first);
jacint@109	299	left.set(w,0);
jacint@109	300	level_list[newlevel]=w;
jacint@109	301	}
jacint@109	302	}
jacint@109	303
jacint@109	304
jacint@109	305	++relabel;
jacint@109	306	if ( relabel >= heur ) {
jacint@109	307	relabel=0;
jacint@109	308	if ( what_heur ) {
jacint@109	309	what_heur=0;
jacint@109	310	heur=heur0;
jacint@109	311	end=false;
jacint@109	312	} else {
jacint@109	313	what_heur=1;
jacint@109	314	heur=heur1;
jacint@109	315	b=k;
jacint@109	316	}
jacint@109	317	}
jacint@109	318
jacint@109	319
jacint@109	320	} // if ( exc > 0 )
jacint@109	321
jacint@109	322
jacint@109	323	} // if stack[b] is nonempty
jacint@109	324
jacint@109	325	} // while(true)
jacint@109	326
jacint@109	327
jacint@109	328
jacint@109	329	for( EachNodeIt v=G.template first<EachNodeIt>();
jacint@109	330	v.valid(); ++v)
jacint@109	331	if (level.get(v) >= n ) cut.set(v,true);
jacint@109	332
jacint@109	333	value = excess.get(t);
jacint@109	334	/Max flow value./
jacint@109	335
jacint@109	336	} //void run()
jacint@109	337
jacint@109	338
jacint@109	339
jacint@109	340
jacint@109	341	T maxFlow() {
jacint@109	342	return value;
jacint@109	343	}
jacint@109	344
jacint@109	345
jacint@109	346
jacint@109	347	CutMap minCut() {
jacint@109	348	return cut;
jacint@109	349	}
jacint@109	350
jacint@109	351
jacint@109	352	};
jacint@109	353	}//namespace
jacint@109	354	#endif
jacint@109	355
jacint@109	356
jacint@109	357
jacint@109	358

author	jacint
	Sun, 22 Feb 2004 12:17:16 +0000
changeset 113	cf7b01232d86
permissions	-rw-r--r--