|
1 /* -*- C++ -*- |
|
2 * lemon/preflow_matching.h - Part of LEMON, a generic C++ optimization library |
|
3 * |
|
4 * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport |
|
5 * (Egervary Research Group on Combinatorial Optimization, EGRES). |
|
6 * |
|
7 * Permission to use, modify and distribute this software is granted |
|
8 * provided that this copyright notice appears in all copies. For |
|
9 * precise terms see the accompanying LICENSE file. |
|
10 * |
|
11 * This software is provided "AS IS" with no warranty of any kind, |
|
12 * express or implied, and with no claim as to its suitability for any |
|
13 * purpose. |
|
14 * |
|
15 */ |
|
16 |
|
17 #ifndef LEMON_CIRCULATION_H |
|
18 #define LEMON_CIRCULATION_H |
|
19 |
|
20 #include <lemon/graph_utils.h> |
|
21 #include <iostream> |
|
22 #include <queue> |
|
23 #include <lemon/tolerance.h> |
|
24 #include <lemon/elevator.h> |
|
25 |
|
26 ///\ingroup flowalgs |
|
27 ///\file |
|
28 ///\brief Push-prelabel algorithm for finding a feasible circulation. |
|
29 /// |
|
30 namespace lemon { |
|
31 |
|
32 ///Preflow algorithm for the Network Circulation Problem. |
|
33 |
|
34 ///\ingroup flowalgs |
|
35 ///This class implements a preflow algorithm |
|
36 ///for the Network Circulation Problem. |
|
37 ///The exact formulation of this problem is the following. |
|
38 /// \f[\sum_{e\in\rho(v)}x(e)-\sum_{e\in\delta(v)}x(e)\leq delta(v)\quad \forall v\in V \f] |
|
39 /// \f[ lo(e)\leq x(e) \leq up(e) \quad \forall e\in E \f] |
|
40 /// |
|
41 template<class Graph, |
|
42 class Value, |
|
43 class FlowMap=typename Graph::template EdgeMap<Value>, |
|
44 class LCapMap=typename Graph::template EdgeMap<Value>, |
|
45 class UCapMap=LCapMap, |
|
46 class DeltaMap=typename Graph::template NodeMap<Value> |
|
47 > |
|
48 class Circulation { |
|
49 typedef typename Graph::Node Node; |
|
50 typedef typename Graph::NodeIt NodeIt; |
|
51 typedef typename Graph::Edge Edge; |
|
52 typedef typename Graph::EdgeIt EdgeIt; |
|
53 typedef typename Graph::InEdgeIt InEdgeIt; |
|
54 typedef typename Graph::OutEdgeIt OutEdgeIt; |
|
55 |
|
56 |
|
57 const Graph &_g; |
|
58 int _node_num; |
|
59 |
|
60 const LCapMap &_lo; |
|
61 const UCapMap &_up; |
|
62 const DeltaMap &_delta; |
|
63 FlowMap &_x; |
|
64 Tolerance<Value> _tol; |
|
65 Elevator<Graph,typename Graph::Node> _levels; |
|
66 typename Graph::template NodeMap<Value> _excess; |
|
67 |
|
68 public: |
|
69 ///\e |
|
70 Circulation(const Graph &g,const LCapMap &lo,const UCapMap &up, |
|
71 const DeltaMap &delta,FlowMap &x) : |
|
72 _g(g), |
|
73 _node_num(countNodes(g)), |
|
74 _lo(lo),_up(up),_delta(delta),_x(x), |
|
75 _levels(g,_node_num), |
|
76 _excess(g) |
|
77 { |
|
78 } |
|
79 |
|
80 private: |
|
81 |
|
82 void addExcess(Node n,Value v) |
|
83 { |
|
84 if(_tol.positive(_excess[n]+=v)) |
|
85 { |
|
86 if(!_levels.active(n)) _levels.activate(n); |
|
87 } |
|
88 else if(_levels.active(n)) _levels.deactivate(n); |
|
89 } |
|
90 |
|
91 void init() |
|
92 { |
|
93 |
|
94 _x=_lo; |
|
95 |
|
96 for(NodeIt n(_g);n!=INVALID;++n) _excess[n]=-_delta[n]; |
|
97 |
|
98 for(EdgeIt e(_g);e!=INVALID;++e) |
|
99 { |
|
100 _excess[_g.target(e)]+=_x[e]; |
|
101 _excess[_g.source(e)]-=_x[e]; |
|
102 } |
|
103 |
|
104 _levels.initStart(); |
|
105 for(NodeIt n(_g);n!=INVALID;++n) |
|
106 _levels.initAddItem(n); |
|
107 _levels.initFinish(); |
|
108 for(NodeIt n(_g);n!=INVALID;++n) |
|
109 if(_tol.positive(_excess[n])) |
|
110 _levels.activate(n); |
|
111 } |
|
112 |
|
113 public: |
|
114 ///Check if \c x is a feasible circulation |
|
115 template<class FT> |
|
116 bool checkX(FT &x) { |
|
117 for(EdgeIt e(_g);e!=INVALID;++e) |
|
118 if(x[e]<_lo[e]||x[e]>_up[e]) return false; |
|
119 for(NodeIt n(_g);n!=INVALID;++n) |
|
120 { |
|
121 Value dif=_delta[n]; |
|
122 for(InEdgeIt e(_g,n);e!=INVALID;++e) dif-=x[e]; |
|
123 for(OutEdgeIt e(_g,n);e!=INVALID;++e) dif+=x[e]; |
|
124 if(_tol.negative(dif)) return false; |
|
125 } |
|
126 return true; |
|
127 }; |
|
128 ///Check if the default \c x is a feasible circulation |
|
129 bool checkX() { return checkX(_x); } |
|
130 |
|
131 ///Chech if \c bar is a real barrier |
|
132 |
|
133 ///Chech if \c bar is a real barrier |
|
134 ///\sa barrier() |
|
135 template<class GT> |
|
136 bool checkBarrier(GT &bar) |
|
137 { |
|
138 Value delta=0; |
|
139 for(NodeIt n(_g);n!=INVALID;++n) |
|
140 if(bar[n]) |
|
141 delta+=_delta[n]; |
|
142 for(EdgeIt e(_g);e!=INVALID;++e) |
|
143 { |
|
144 Node s=_g.source(e); |
|
145 Node t=_g.target(e); |
|
146 if(bar[s]&&!bar[t]) delta+=_up[e]; |
|
147 else if(bar[t]&&!bar[s]) delta-=_lo[e]; |
|
148 } |
|
149 return _tol.negative(delta); |
|
150 } |
|
151 ///Chech whether or not the last execution provides a barrier |
|
152 |
|
153 ///Chech whether or not the last execution provides a barrier |
|
154 ///\sa barrier() |
|
155 bool checkBarrier() |
|
156 { |
|
157 typename Graph:: template NodeMap<bool> bar(_g); |
|
158 barrier(bar); |
|
159 return checkBarrier(bar); |
|
160 } |
|
161 ///Run the algorithm |
|
162 |
|
163 ///This function runs the algorithm. |
|
164 ///\return This function returns -1 if it found a feasible circulation. |
|
165 /// nonnegative values (including 0) mean that no feasible solution is |
|
166 /// found. In this case the return value means an "empty level". |
|
167 /// |
|
168 ///\sa barrier() |
|
169 int run() |
|
170 { |
|
171 init(); |
|
172 |
|
173 #ifdef LEMON_CIRCULATION_DEBUG |
|
174 for(NodeIt n(_g);n!=INVALID;++n) |
|
175 std::cerr<< _levels[n] << ' '; |
|
176 std::cerr << std::endl; |
|
177 #endif |
|
178 Node act; |
|
179 Node bact=INVALID; |
|
180 Node last_activated=INVALID; |
|
181 while((act=_levels.highestActive())!=INVALID) { |
|
182 int actlevel=_levels[act]; |
|
183 int tlevel; |
|
184 int mlevel=_node_num; |
|
185 Value exc=_excess[act]; |
|
186 Value fc; |
|
187 |
|
188 #ifdef LEMON_CIRCULATION_DEBUG |
|
189 for(NodeIt n(_g);n!=INVALID;++n) |
|
190 std::cerr<< _levels[n] << ' '; |
|
191 std::cerr << std::endl; |
|
192 std::cerr << "Process node " << _g.id(act) |
|
193 << " on level " << actlevel |
|
194 << " with excess " << exc |
|
195 << std::endl; |
|
196 #endif |
|
197 for(OutEdgeIt e(_g,act);e!=INVALID; ++e) |
|
198 if((fc=_up[e]-_x[e])>0) |
|
199 if((tlevel=_levels[_g.target(e)])<actlevel) |
|
200 if(fc<=exc) { |
|
201 _x[e]=_up[e]; |
|
202 addExcess(_g.target(e),fc); |
|
203 exc-=fc; |
|
204 #ifdef LEMON_CIRCULATION_DEBUG |
|
205 std::cerr << " Push " << fc |
|
206 << " toward " << _g.id(_g.target(e)) << std::endl; |
|
207 #endif |
|
208 } |
|
209 else { |
|
210 _x[e]+=exc; |
|
211 addExcess(_g.target(e),exc); |
|
212 //exc=0; |
|
213 _excess[act]=0; |
|
214 _levels.deactivate(act); |
|
215 #ifdef LEMON_CIRCULATION_DEBUG |
|
216 std::cerr << " Push " << exc |
|
217 << " toward " << _g.id(_g.target(e)) << std::endl; |
|
218 std::cerr << " Deactivate\n"; |
|
219 #endif |
|
220 goto next_l; |
|
221 } |
|
222 else if(tlevel<mlevel) mlevel=tlevel; |
|
223 |
|
224 for(InEdgeIt e(_g,act);e!=INVALID; ++e) |
|
225 if((fc=_x[e]-_lo[e])>0) |
|
226 if((tlevel=_levels[_g.source(e)])<actlevel) |
|
227 if(fc<=exc) { |
|
228 _x[e]=_lo[e]; |
|
229 addExcess(_g.source(e),fc); |
|
230 exc-=fc; |
|
231 #ifdef LEMON_CIRCULATION_DEBUG |
|
232 std::cerr << " Push " << fc |
|
233 << " toward " << _g.id(_g.source(e)) << std::endl; |
|
234 #endif |
|
235 } |
|
236 else { |
|
237 _x[e]-=exc; |
|
238 addExcess(_g.source(e),exc); |
|
239 //exc=0; |
|
240 _excess[act]=0; |
|
241 _levels.deactivate(act); |
|
242 #ifdef LEMON_CIRCULATION_DEBUG |
|
243 std::cerr << " Push " << exc |
|
244 << " toward " << _g.id(_g.source(e)) << std::endl; |
|
245 std::cerr << " Deactivate\n"; |
|
246 #endif |
|
247 goto next_l; |
|
248 } |
|
249 else if(tlevel<mlevel) mlevel=tlevel; |
|
250 |
|
251 _excess[act]=exc; |
|
252 if(!_tol.positive(exc)) _levels.deactivate(act); |
|
253 else if(mlevel==_node_num) { |
|
254 _levels.liftHighestActiveTo(_node_num); |
|
255 #ifdef LEMON_CIRCULATION_DEBUG |
|
256 std::cerr << " Lift to level " << _node_num << std::endl; |
|
257 #endif |
|
258 return _levels.onLevel(_node_num-1)==0?_node_num-1:actlevel; |
|
259 } |
|
260 else { |
|
261 _levels.liftHighestActiveTo(mlevel+1); |
|
262 #ifdef LEMON_CIRCULATION_DEBUG |
|
263 std::cerr << " Lift to level " << mlevel+1 << std::endl; |
|
264 #endif |
|
265 if(_levels.onLevel(actlevel)==0) |
|
266 return actlevel; |
|
267 } |
|
268 next_l: |
|
269 ; |
|
270 } |
|
271 #ifdef LEMON_CIRCULATION_DEBUG |
|
272 std::cerr << "Feasible flow found.\n"; |
|
273 #endif |
|
274 return -1; |
|
275 } |
|
276 |
|
277 ///Return a barrier |
|
278 |
|
279 ///Barrier is a set \e B of nodes for which |
|
280 /// \f[ \sum_{v\in B}delta(v)<\sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f] |
|
281 ///holds. The existence of a set with this property prooves that a feasible |
|
282 ///flow cannot exists. |
|
283 ///\pre The run() must have been executed, and its return value was -1. |
|
284 ///\sa checkBarrier() |
|
285 ///\sa run() |
|
286 template<class GT> |
|
287 void barrier(GT &bar,int empty_level=-1) |
|
288 { |
|
289 if(empty_level==-1) |
|
290 for(empty_level=0;_levels.onLevel(empty_level);empty_level++) ; |
|
291 for(NodeIt n(_g);n!=INVALID;++n) |
|
292 bar[n] = _levels[n]>empty_level; |
|
293 } |
|
294 |
|
295 }; |
|
296 |
|
297 } |
|
298 |
|
299 #endif |