154 /// \brief Constructor |
127 /// \brief Constructor |
155 /// |
128 /// |
156 /// Constructor of the algorithm class. |
129 /// Constructor of the algorithm class. |
157 HaoOrlin(const Graph& graph, const CapacityMap& capacity, |
130 HaoOrlin(const Graph& graph, const CapacityMap& capacity, |
158 const Tolerance& tolerance = Tolerance()) : |
131 const Tolerance& tolerance = Tolerance()) : |
159 _graph(&graph), _capacity(&capacity), |
132 _graph(graph), _capacity(&capacity), _flow(0), _source(), |
160 _preflow(0), _source(), _target(), |
133 _node_num(), _first(), _last(), _next(0), _prev(0), |
161 _out_res_graph(0), _rev_graph(0), _in_res_graph(0), |
134 _active(0), _bucket(0), _dormant(), _sets(), _highest(), |
162 _wake(0),_dist(0), _excess(0), _source_set(0), |
135 _excess(0), _source_set(0), _min_cut(), _min_cut_map(0), |
163 _highest_active(), _active_nodes(), _dormant_max(), _dormant(), |
136 _tolerance(tolerance) {} |
164 _min_cut(), _min_cut_map(0), _tolerance(tolerance) {} |
|
165 |
137 |
166 ~HaoOrlin() { |
138 ~HaoOrlin() { |
167 if (_min_cut_map) { |
139 if (_min_cut_map) { |
168 delete _min_cut_map; |
140 delete _min_cut_map; |
169 } |
141 } |
170 if (_in_res_graph) { |
|
171 delete _in_res_graph; |
|
172 } |
|
173 if (_rev_graph) { |
|
174 delete _rev_graph; |
|
175 } |
|
176 if (_out_res_graph) { |
|
177 delete _out_res_graph; |
|
178 } |
|
179 if (_source_set) { |
142 if (_source_set) { |
180 delete _source_set; |
143 delete _source_set; |
181 } |
144 } |
182 if (_excess) { |
145 if (_excess) { |
183 delete _excess; |
146 delete _excess; |
184 } |
147 } |
185 if (_dist) { |
148 if (_next) { |
186 delete _dist; |
149 delete _next; |
187 } |
150 } |
188 if (_wake) { |
151 if (_prev) { |
189 delete _wake; |
152 delete _prev; |
190 } |
153 } |
191 if (_preflow) { |
154 if (_active) { |
192 delete _preflow; |
155 delete _active; |
|
156 } |
|
157 if (_bucket) { |
|
158 delete _bucket; |
|
159 } |
|
160 if (_flow) { |
|
161 delete _flow; |
193 } |
162 } |
194 } |
163 } |
195 |
164 |
196 private: |
165 private: |
|
166 |
|
167 void activate(const Node& i) { |
|
168 _active->set(i, true); |
|
169 |
|
170 int bucket = (*_bucket)[i]; |
|
171 |
|
172 if ((*_prev)[i] == INVALID || (*_active)[(*_prev)[i]]) return; |
|
173 //unlace |
|
174 _next->set((*_prev)[i], (*_next)[i]); |
|
175 if ((*_next)[i] != INVALID) { |
|
176 _prev->set((*_next)[i], (*_prev)[i]); |
|
177 } else { |
|
178 _last[bucket] = (*_prev)[i]; |
|
179 } |
|
180 //lace |
|
181 _next->set(i, _first[bucket]); |
|
182 _prev->set(_first[bucket], i); |
|
183 _prev->set(i, INVALID); |
|
184 _first[bucket] = i; |
|
185 } |
|
186 |
|
187 void deactivate(const Node& i) { |
|
188 _active->set(i, false); |
|
189 int bucket = (*_bucket)[i]; |
|
190 |
|
191 if ((*_next)[i] == INVALID || !(*_active)[(*_next)[i]]) return; |
|
192 |
|
193 //unlace |
|
194 _prev->set((*_next)[i], (*_prev)[i]); |
|
195 if ((*_prev)[i] != INVALID) { |
|
196 _next->set((*_prev)[i], (*_next)[i]); |
|
197 } else { |
|
198 _first[bucket] = (*_next)[i]; |
|
199 } |
|
200 //lace |
|
201 _prev->set(i, _last[bucket]); |
|
202 _next->set(_last[bucket], i); |
|
203 _next->set(i, INVALID); |
|
204 _last[bucket] = i; |
|
205 } |
|
206 |
|
207 void addItem(const Node& i, int bucket) { |
|
208 (*_bucket)[i] = bucket; |
|
209 if (_last[bucket] != INVALID) { |
|
210 _prev->set(i, _last[bucket]); |
|
211 _next->set(_last[bucket], i); |
|
212 _next->set(i, INVALID); |
|
213 _last[bucket] = i; |
|
214 } else { |
|
215 _prev->set(i, INVALID); |
|
216 _first[bucket] = i; |
|
217 _next->set(i, INVALID); |
|
218 _last[bucket] = i; |
|
219 } |
|
220 } |
197 |
221 |
198 template <typename ResGraph> |
222 void findMinCutOut() { |
199 void findMinCut(const Node& target, bool out, ResGraph& res_graph) { |
223 |
200 typedef typename ResGraph::Edge ResEdge; |
224 for (NodeIt n(_graph); n != INVALID; ++n) { |
201 typedef typename ResGraph::OutEdgeIt ResOutEdgeIt; |
225 _excess->set(n, 0); |
202 |
226 } |
203 for (typename Graph::EdgeIt it(*_graph); it != INVALID; ++it) { |
227 |
204 (*_preflow)[it] = 0; |
228 for (EdgeIt e(_graph); e != INVALID; ++e) { |
205 } |
229 _flow->set(e, 0); |
206 for (NodeIt it(*_graph); it != INVALID; ++it) { |
230 } |
207 (*_wake)[it] = true; |
231 |
208 (*_dist)[it] = 1; |
232 int bucket_num = 1; |
209 (*_excess)[it] = 0; |
233 |
210 (*_source_set)[it] = false; |
234 { |
211 } |
235 typename Graph::template NodeMap<bool> reached(_graph, false); |
212 |
236 |
213 _dormant[0].push_front(_source); |
237 reached.set(_source, true); |
214 (*_source_set)[_source] = true; |
238 |
215 _dormant_max = 0; |
239 bool first_set = true; |
216 (*_wake)[_source] = false; |
240 |
217 |
241 for (NodeIt t(_graph); t != INVALID; ++t) { |
218 _level_size[0] = 1; |
242 if (reached[t]) continue; |
219 _level_size[1] = _node_num - 1; |
243 _sets.push_front(std::list<int>()); |
220 |
244 _sets.front().push_front(bucket_num); |
221 _target = target; |
245 _dormant[bucket_num] = !first_set; |
222 (*_dist)[target] = 0; |
246 |
223 |
247 _bucket->set(t, bucket_num); |
224 for (ResOutEdgeIt it(res_graph, _source); it != INVALID; ++it) { |
248 _first[bucket_num] = _last[bucket_num] = t; |
225 Value delta = res_graph.rescap(it); |
249 _next->set(t, INVALID); |
226 (*_excess)[_source] -= delta; |
250 _prev->set(t, INVALID); |
227 res_graph.augment(it, delta); |
251 |
228 Node a = res_graph.target(it); |
252 ++bucket_num; |
229 if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) { |
253 |
230 _active_nodes[(*_dist)[a]].push_front(a); |
254 std::vector<Node> queue; |
231 if (_highest_active < (*_dist)[a]) { |
255 queue.push_back(t); |
232 _highest_active = (*_dist)[a]; |
256 reached.set(t, true); |
233 } |
257 |
234 } |
258 while (!queue.empty()) { |
235 (*_excess)[a] += delta; |
259 _sets.front().push_front(bucket_num); |
236 } |
260 _dormant[bucket_num] = !first_set; |
237 |
261 _first[bucket_num] = _last[bucket_num] = INVALID; |
238 |
262 |
239 do { |
263 std::vector<Node> nqueue; |
240 Node n; |
264 for (int i = 0; i < int(queue.size()); ++i) { |
241 while ((n = findActiveNode()) != INVALID) { |
265 Node n = queue[i]; |
242 for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) { |
266 for (InEdgeIt e(_graph, n); e != INVALID; ++e) { |
243 Node a = res_graph.target(e); |
267 Node u = _graph.source(e); |
244 if ((*_dist)[a] >= (*_dist)[n] || !(*_wake)[a]) continue; |
268 if (!reached[u] && _tolerance.positive((*_capacity)[e])) { |
245 Value delta = res_graph.rescap(e); |
269 reached.set(u, true); |
246 if (_tolerance.positive((*_excess)[n] - delta)) { |
270 addItem(u, bucket_num); |
247 (*_excess)[n] -= delta; |
271 nqueue.push_back(u); |
|
272 } |
|
273 } |
|
274 } |
|
275 queue.swap(nqueue); |
|
276 ++bucket_num; |
|
277 } |
|
278 _sets.front().pop_front(); |
|
279 --bucket_num; |
|
280 first_set = false; |
|
281 } |
|
282 |
|
283 _bucket->set(_source, 0); |
|
284 _dormant[0] = true; |
|
285 } |
|
286 _source_set->set(_source, true); |
|
287 |
|
288 Node target = _last[_sets.back().back()]; |
|
289 { |
|
290 for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) { |
|
291 if (_tolerance.positive((*_capacity)[e])) { |
|
292 Node u = _graph.target(e); |
|
293 _flow->set(e, (*_capacity)[e]); |
|
294 _excess->set(u, (*_excess)[u] + (*_capacity)[e]); |
|
295 if (!(*_active)[u] && u != _source) { |
|
296 activate(u); |
|
297 } |
|
298 } |
|
299 } |
|
300 if ((*_active)[target]) { |
|
301 deactivate(target); |
|
302 } |
|
303 |
|
304 _highest = _sets.back().begin(); |
|
305 while (_highest != _sets.back().end() && |
|
306 !(*_active)[_first[*_highest]]) { |
|
307 ++_highest; |
|
308 } |
|
309 } |
|
310 |
|
311 |
|
312 while (true) { |
|
313 while (_highest != _sets.back().end()) { |
|
314 Node n = _first[*_highest]; |
|
315 Value excess = (*_excess)[n]; |
|
316 int next_bucket = _node_num; |
|
317 |
|
318 int under_bucket; |
|
319 if (++std::list<int>::iterator(_highest) == _sets.back().end()) { |
|
320 under_bucket = -1; |
|
321 } else { |
|
322 under_bucket = *(++std::list<int>::iterator(_highest)); |
|
323 } |
|
324 |
|
325 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) { |
|
326 Node v = _graph.target(e); |
|
327 if (_dormant[(*_bucket)[v]]) continue; |
|
328 Value rem = (*_capacity)[e] - (*_flow)[e]; |
|
329 if (!_tolerance.positive(rem)) continue; |
|
330 if ((*_bucket)[v] == under_bucket) { |
|
331 if (!(*_active)[v] && v != target) { |
|
332 activate(v); |
|
333 } |
|
334 if (!_tolerance.less(rem, excess)) { |
|
335 _flow->set(e, (*_flow)[e] + excess); |
|
336 _excess->set(v, (*_excess)[v] + excess); |
|
337 excess = 0; |
|
338 goto no_more_push; |
|
339 } else { |
|
340 excess -= rem; |
|
341 _excess->set(v, (*_excess)[v] + rem); |
|
342 _flow->set(e, (*_capacity)[e]); |
|
343 } |
|
344 } else if (next_bucket > (*_bucket)[v]) { |
|
345 next_bucket = (*_bucket)[v]; |
|
346 } |
|
347 } |
|
348 |
|
349 for (InEdgeIt e(_graph, n); e != INVALID; ++e) { |
|
350 Node v = _graph.source(e); |
|
351 if (_dormant[(*_bucket)[v]]) continue; |
|
352 Value rem = (*_flow)[e]; |
|
353 if (!_tolerance.positive(rem)) continue; |
|
354 if ((*_bucket)[v] == under_bucket) { |
|
355 if (!(*_active)[v] && v != target) { |
|
356 activate(v); |
|
357 } |
|
358 if (!_tolerance.less(rem, excess)) { |
|
359 _flow->set(e, (*_flow)[e] - excess); |
|
360 _excess->set(v, (*_excess)[v] + excess); |
|
361 excess = 0; |
|
362 goto no_more_push; |
|
363 } else { |
|
364 excess -= rem; |
|
365 _excess->set(v, (*_excess)[v] + rem); |
|
366 _flow->set(e, 0); |
|
367 } |
|
368 } else if (next_bucket > (*_bucket)[v]) { |
|
369 next_bucket = (*_bucket)[v]; |
|
370 } |
|
371 } |
|
372 |
|
373 no_more_push: |
|
374 |
|
375 _excess->set(n, excess); |
|
376 |
|
377 if (excess != 0) { |
|
378 if ((*_next)[n] == INVALID) { |
|
379 typename std::list<std::list<int> >::iterator new_set = |
|
380 _sets.insert(--_sets.end(), std::list<int>()); |
|
381 new_set->splice(new_set->end(), _sets.back(), |
|
382 _sets.back().begin(), ++_highest); |
|
383 for (std::list<int>::iterator it = new_set->begin(); |
|
384 it != new_set->end(); ++it) { |
|
385 _dormant[*it] = true; |
|
386 } |
|
387 while (_highest != _sets.back().end() && |
|
388 !(*_active)[_first[*_highest]]) { |
|
389 ++_highest; |
|
390 } |
|
391 } else if (next_bucket == _node_num) { |
|
392 _first[(*_bucket)[n]] = (*_next)[n]; |
|
393 _prev->set((*_next)[n], INVALID); |
|
394 |
|
395 std::list<std::list<int> >::iterator new_set = |
|
396 _sets.insert(--_sets.end(), std::list<int>()); |
|
397 |
|
398 new_set->push_front(bucket_num); |
|
399 _bucket->set(n, bucket_num); |
|
400 _first[bucket_num] = _last[bucket_num] = n; |
|
401 _next->set(n, INVALID); |
|
402 _prev->set(n, INVALID); |
|
403 _dormant[bucket_num] = true; |
|
404 ++bucket_num; |
|
405 |
|
406 while (_highest != _sets.back().end() && |
|
407 !(*_active)[_first[*_highest]]) { |
|
408 ++_highest; |
|
409 } |
248 } else { |
410 } else { |
249 delta = (*_excess)[n]; |
411 _first[*_highest] = (*_next)[n]; |
250 (*_excess)[n] = 0; |
412 _prev->set((*_next)[n], INVALID); |
251 } |
413 |
252 res_graph.augment(e, delta); |
414 while (next_bucket != *_highest) { |
253 if ((*_excess)[a] == 0 && a != _target) { |
415 --_highest; |
254 _active_nodes[(*_dist)[a]].push_front(a); |
416 } |
255 } |
417 |
256 (*_excess)[a] += delta; |
418 if (_highest == _sets.back().begin()) { |
257 if ((*_excess)[n] == 0) break; |
419 _sets.back().push_front(bucket_num); |
258 } |
420 _dormant[bucket_num] = false; |
259 if ((*_excess)[n] != 0) { |
421 _first[bucket_num] = _last[bucket_num] = INVALID; |
260 relabel(n, res_graph); |
422 ++bucket_num; |
261 } |
423 } |
262 } |
424 --_highest; |
263 |
425 |
264 Value current_value = cutValue(out); |
426 _bucket->set(n, *_highest); |
265 if (_min_cut > current_value){ |
427 _next->set(n, _first[*_highest]); |
266 if (out) { |
428 if (_first[*_highest] != INVALID) { |
267 for (NodeIt it(*_graph); it != INVALID; ++it) { |
429 _prev->set(_first[*_highest], n); |
268 _min_cut_map->set(it, !(*_wake)[it]); |
430 } else { |
269 } |
431 _last[*_highest] = n; |
270 } else { |
432 } |
271 for (NodeIt it(*_graph); it != INVALID; ++it) { |
433 _first[*_highest] = n; |
272 _min_cut_map->set(it, (*_wake)[it]); |
434 } |
273 } |
435 } else { |
274 } |
436 |
275 |
437 deactivate(n); |
276 _min_cut = current_value; |
438 if (!(*_active)[_first[*_highest]]) { |
277 } |
439 ++_highest; |
278 |
440 if (_highest != _sets.back().end() && |
279 } while (selectNewSink(res_graph)); |
441 !(*_active)[_first[*_highest]]) { |
280 } |
442 _highest = _sets.back().end(); |
281 |
443 } |
282 template <typename ResGraph> |
444 } |
283 void relabel(const Node& n, ResGraph& res_graph) { |
445 } |
284 typedef typename ResGraph::OutEdgeIt ResOutEdgeIt; |
446 } |
285 |
447 |
286 int k = (*_dist)[n]; |
448 if ((*_excess)[target] < _min_cut) { |
287 if (_level_size[k] == 1) { |
449 _min_cut = (*_excess)[target]; |
288 ++_dormant_max; |
450 for (NodeIt i(_graph); i != INVALID; ++i) { |
289 for (NodeIt it(*_graph); it != INVALID; ++it) { |
451 _min_cut_map->set(i, true); |
290 if ((*_wake)[it] && (*_dist)[it] >= k) { |
452 } |
291 (*_wake)[it] = false; |
453 for (std::list<int>::iterator it = _sets.back().begin(); |
292 _dormant[_dormant_max].push_front(it); |
454 it != _sets.back().end(); ++it) { |
293 --_level_size[(*_dist)[it]]; |
455 Node n = _first[*it]; |
294 } |
456 while (n != INVALID) { |
295 } |
457 _min_cut_map->set(n, false); |
296 --_highest_active; |
458 n = (*_next)[n]; |
297 } else { |
459 } |
298 int new_dist = _node_num; |
460 } |
299 for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) { |
461 } |
300 Node t = res_graph.target(e); |
462 |
301 if ((*_wake)[t] && new_dist > (*_dist)[t]) { |
463 { |
302 new_dist = (*_dist)[t]; |
464 Node new_target; |
303 } |
465 if ((*_prev)[target] != INVALID) { |
304 } |
466 _last[(*_bucket)[target]] = (*_prev)[target]; |
305 if (new_dist == _node_num) { |
467 _next->set((*_prev)[target], INVALID); |
306 ++_dormant_max; |
468 new_target = (*_prev)[target]; |
307 (*_wake)[n] = false; |
469 } else { |
308 _dormant[_dormant_max].push_front(n); |
470 _sets.back().pop_back(); |
309 --_level_size[(*_dist)[n]]; |
471 if (_sets.back().empty()) { |
310 } else { |
472 _sets.pop_back(); |
311 --_level_size[(*_dist)[n]]; |
473 if (_sets.empty()) |
312 (*_dist)[n] = new_dist + 1; |
474 break; |
313 _highest_active = (*_dist)[n]; |
475 for (std::list<int>::iterator it = _sets.back().begin(); |
314 _active_nodes[_highest_active].push_front(n); |
476 it != _sets.back().end(); ++it) { |
315 ++_level_size[(*_dist)[n]]; |
477 _dormant[*it] = false; |
316 } |
478 } |
317 } |
479 } |
318 } |
480 new_target = _last[_sets.back().back()]; |
319 |
481 } |
320 template <typename ResGraph> |
482 |
321 bool selectNewSink(ResGraph& res_graph) { |
483 _bucket->set(target, 0); |
322 typedef typename ResGraph::OutEdgeIt ResOutEdgeIt; |
484 |
323 |
485 _source_set->set(target, true); |
324 Node old_target = _target; |
486 for (OutEdgeIt e(_graph, target); e != INVALID; ++e) { |
325 (*_wake)[_target] = false; |
487 Value rem = (*_capacity)[e] - (*_flow)[e]; |
326 --_level_size[(*_dist)[_target]]; |
488 if (!_tolerance.positive(rem)) continue; |
327 _dormant[0].push_front(_target); |
489 Node v = _graph.target(e); |
328 (*_source_set)[_target] = true; |
490 if (!(*_active)[v] && !(*_source_set)[v]) { |
329 if (int(_dormant[0].size()) == _node_num){ |
491 activate(v); |
330 _dormant[0].clear(); |
492 } |
331 return false; |
493 _excess->set(v, (*_excess)[v] + rem); |
332 } |
494 _flow->set(e, (*_capacity)[e]); |
333 |
495 } |
334 bool wake_was_empty = false; |
496 |
335 |
497 for (InEdgeIt e(_graph, target); e != INVALID; ++e) { |
336 if(_wake->trueNum() == 0) { |
498 Value rem = (*_flow)[e]; |
337 while (!_dormant[_dormant_max].empty()){ |
499 if (!_tolerance.positive(rem)) continue; |
338 (*_wake)[_dormant[_dormant_max].front()] = true; |
500 Node v = _graph.source(e); |
339 ++_level_size[(*_dist)[_dormant[_dormant_max].front()]]; |
501 if (!(*_active)[v] && !(*_source_set)[v]) { |
340 _dormant[_dormant_max].pop_front(); |
502 activate(v); |
341 } |
503 } |
342 --_dormant_max; |
504 _excess->set(v, (*_excess)[v] + rem); |
343 wake_was_empty = true; |
505 _flow->set(e, 0); |
344 } |
506 } |
345 |
507 |
346 int min_dist = std::numeric_limits<int>::max(); |
508 target = new_target; |
347 for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) { |
509 if ((*_active)[target]) { |
348 if (min_dist > (*_dist)[it]){ |
510 deactivate(target); |
349 _target = it; |
511 } |
350 min_dist = (*_dist)[it]; |
512 |
351 } |
513 _highest = _sets.back().begin(); |
352 } |
514 while (_highest != _sets.back().end() && |
353 |
515 !(*_active)[_first[*_highest]]) { |
354 if (wake_was_empty){ |
516 ++_highest; |
355 for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) { |
517 } |
356 if ((*_excess)[it] != 0 && it != _target) { |
518 } |
357 _active_nodes[(*_dist)[it]].push_front(it); |
519 } |
358 if (_highest_active < (*_dist)[it]) { |
520 } |
359 _highest_active = (*_dist)[it]; |
521 |
360 } |
522 void findMinCutIn() { |
361 } |
523 |
362 } |
524 for (NodeIt n(_graph); n != INVALID; ++n) { |
363 } |
525 _excess->set(n, 0); |
364 |
526 } |
365 Node n = old_target; |
527 |
366 for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e){ |
528 for (EdgeIt e(_graph); e != INVALID; ++e) { |
367 Node a = res_graph.target(e); |
529 _flow->set(e, 0); |
368 if (!(*_wake)[a]) continue; |
530 } |
369 Value delta = res_graph.rescap(e); |
531 |
370 res_graph.augment(e, delta); |
532 int bucket_num = 1; |
371 (*_excess)[n] -= delta; |
|
372 if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) { |
|
373 _active_nodes[(*_dist)[a]].push_front(a); |
|
374 if (_highest_active < (*_dist)[a]) { |
|
375 _highest_active = (*_dist)[a]; |
|
376 } |
|
377 } |
|
378 (*_excess)[a] += delta; |
|
379 } |
|
380 |
533 |
381 return true; |
534 { |
382 } |
535 typename Graph::template NodeMap<bool> reached(_graph, false); |
383 |
536 |
384 Node findActiveNode() { |
537 reached.set(_source, true); |
385 while (_highest_active > 0 && _active_nodes[_highest_active].empty()){ |
538 |
386 --_highest_active; |
539 bool first_set = true; |
387 } |
540 |
388 if( _highest_active > 0) { |
541 for (NodeIt t(_graph); t != INVALID; ++t) { |
389 Node n = _active_nodes[_highest_active].front(); |
542 if (reached[t]) continue; |
390 _active_nodes[_highest_active].pop_front(); |
543 _sets.push_front(std::list<int>()); |
391 return n; |
544 _sets.front().push_front(bucket_num); |
392 } else { |
545 _dormant[bucket_num] = !first_set; |
393 return INVALID; |
546 |
394 } |
547 _bucket->set(t, bucket_num); |
395 } |
548 _first[bucket_num] = _last[bucket_num] = t; |
396 |
549 _next->set(t, INVALID); |
397 Value cutValue(bool out) { |
550 _prev->set(t, INVALID); |
398 Value value = 0; |
551 |
399 if (out) { |
552 ++bucket_num; |
400 for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) { |
553 |
401 for (InEdgeIt e(*_graph, it); e != INVALID; ++e) { |
554 std::vector<Node> queue; |
402 if (!(*_wake)[_graph->source(e)]){ |
555 queue.push_back(t); |
403 value += (*_capacity)[e]; |
556 reached.set(t, true); |
404 } |
557 |
405 } |
558 while (!queue.empty()) { |
406 } |
559 _sets.front().push_front(bucket_num); |
407 } else { |
560 _dormant[bucket_num] = !first_set; |
408 for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) { |
561 _first[bucket_num] = _last[bucket_num] = INVALID; |
409 for (OutEdgeIt e(*_graph, it); e != INVALID; ++e) { |
562 |
410 if (!(*_wake)[_graph->target(e)]){ |
563 std::vector<Node> nqueue; |
411 value += (*_capacity)[e]; |
564 for (int i = 0; i < int(queue.size()); ++i) { |
412 } |
565 Node n = queue[i]; |
413 } |
566 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) { |
414 } |
567 Node u = _graph.target(e); |
415 } |
568 if (!reached[u] && _tolerance.positive((*_capacity)[e])) { |
416 return value; |
569 reached.set(u, true); |
417 } |
570 addItem(u, bucket_num); |
418 |
571 nqueue.push_back(u); |
|
572 } |
|
573 } |
|
574 } |
|
575 queue.swap(nqueue); |
|
576 ++bucket_num; |
|
577 } |
|
578 _sets.front().pop_front(); |
|
579 --bucket_num; |
|
580 first_set = false; |
|
581 } |
|
582 |
|
583 _bucket->set(_source, 0); |
|
584 _dormant[0] = true; |
|
585 } |
|
586 _source_set->set(_source, true); |
|
587 |
|
588 Node target = _last[_sets.back().back()]; |
|
589 { |
|
590 for (InEdgeIt e(_graph, _source); e != INVALID; ++e) { |
|
591 if (_tolerance.positive((*_capacity)[e])) { |
|
592 Node u = _graph.source(e); |
|
593 _flow->set(e, (*_capacity)[e]); |
|
594 _excess->set(u, (*_excess)[u] + (*_capacity)[e]); |
|
595 if (!(*_active)[u] && u != _source) { |
|
596 activate(u); |
|
597 } |
|
598 } |
|
599 } |
|
600 if ((*_active)[target]) { |
|
601 deactivate(target); |
|
602 } |
|
603 |
|
604 _highest = _sets.back().begin(); |
|
605 while (_highest != _sets.back().end() && |
|
606 !(*_active)[_first[*_highest]]) { |
|
607 ++_highest; |
|
608 } |
|
609 } |
|
610 |
|
611 |
|
612 while (true) { |
|
613 while (_highest != _sets.back().end()) { |
|
614 Node n = _first[*_highest]; |
|
615 Value excess = (*_excess)[n]; |
|
616 int next_bucket = _node_num; |
|
617 |
|
618 int under_bucket; |
|
619 if (++std::list<int>::iterator(_highest) == _sets.back().end()) { |
|
620 under_bucket = -1; |
|
621 } else { |
|
622 under_bucket = *(++std::list<int>::iterator(_highest)); |
|
623 } |
|
624 |
|
625 for (InEdgeIt e(_graph, n); e != INVALID; ++e) { |
|
626 Node v = _graph.source(e); |
|
627 if (_dormant[(*_bucket)[v]]) continue; |
|
628 Value rem = (*_capacity)[e] - (*_flow)[e]; |
|
629 if (!_tolerance.positive(rem)) continue; |
|
630 if ((*_bucket)[v] == under_bucket) { |
|
631 if (!(*_active)[v] && v != target) { |
|
632 activate(v); |
|
633 } |
|
634 if (!_tolerance.less(rem, excess)) { |
|
635 _flow->set(e, (*_flow)[e] + excess); |
|
636 _excess->set(v, (*_excess)[v] + excess); |
|
637 excess = 0; |
|
638 goto no_more_push; |
|
639 } else { |
|
640 excess -= rem; |
|
641 _excess->set(v, (*_excess)[v] + rem); |
|
642 _flow->set(e, (*_capacity)[e]); |
|
643 } |
|
644 } else if (next_bucket > (*_bucket)[v]) { |
|
645 next_bucket = (*_bucket)[v]; |
|
646 } |
|
647 } |
|
648 |
|
649 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) { |
|
650 Node v = _graph.target(e); |
|
651 if (_dormant[(*_bucket)[v]]) continue; |
|
652 Value rem = (*_flow)[e]; |
|
653 if (!_tolerance.positive(rem)) continue; |
|
654 if ((*_bucket)[v] == under_bucket) { |
|
655 if (!(*_active)[v] && v != target) { |
|
656 activate(v); |
|
657 } |
|
658 if (!_tolerance.less(rem, excess)) { |
|
659 _flow->set(e, (*_flow)[e] - excess); |
|
660 _excess->set(v, (*_excess)[v] + excess); |
|
661 excess = 0; |
|
662 goto no_more_push; |
|
663 } else { |
|
664 excess -= rem; |
|
665 _excess->set(v, (*_excess)[v] + rem); |
|
666 _flow->set(e, 0); |
|
667 } |
|
668 } else if (next_bucket > (*_bucket)[v]) { |
|
669 next_bucket = (*_bucket)[v]; |
|
670 } |
|
671 } |
|
672 |
|
673 no_more_push: |
|
674 |
|
675 _excess->set(n, excess); |
|
676 |
|
677 if (excess != 0) { |
|
678 if ((*_next)[n] == INVALID) { |
|
679 typename std::list<std::list<int> >::iterator new_set = |
|
680 _sets.insert(--_sets.end(), std::list<int>()); |
|
681 new_set->splice(new_set->end(), _sets.back(), |
|
682 _sets.back().begin(), ++_highest); |
|
683 for (std::list<int>::iterator it = new_set->begin(); |
|
684 it != new_set->end(); ++it) { |
|
685 _dormant[*it] = true; |
|
686 } |
|
687 while (_highest != _sets.back().end() && |
|
688 !(*_active)[_first[*_highest]]) { |
|
689 ++_highest; |
|
690 } |
|
691 } else if (next_bucket == _node_num) { |
|
692 _first[(*_bucket)[n]] = (*_next)[n]; |
|
693 _prev->set((*_next)[n], INVALID); |
|
694 |
|
695 std::list<std::list<int> >::iterator new_set = |
|
696 _sets.insert(--_sets.end(), std::list<int>()); |
|
697 |
|
698 new_set->push_front(bucket_num); |
|
699 _bucket->set(n, bucket_num); |
|
700 _first[bucket_num] = _last[bucket_num] = n; |
|
701 _next->set(n, INVALID); |
|
702 _prev->set(n, INVALID); |
|
703 _dormant[bucket_num] = true; |
|
704 ++bucket_num; |
|
705 |
|
706 while (_highest != _sets.back().end() && |
|
707 !(*_active)[_first[*_highest]]) { |
|
708 ++_highest; |
|
709 } |
|
710 } else { |
|
711 _first[*_highest] = (*_next)[n]; |
|
712 _prev->set((*_next)[n], INVALID); |
|
713 |
|
714 while (next_bucket != *_highest) { |
|
715 --_highest; |
|
716 } |
|
717 if (_highest == _sets.back().begin()) { |
|
718 _sets.back().push_front(bucket_num); |
|
719 _dormant[bucket_num] = false; |
|
720 _first[bucket_num] = _last[bucket_num] = INVALID; |
|
721 ++bucket_num; |
|
722 } |
|
723 --_highest; |
|
724 |
|
725 _bucket->set(n, *_highest); |
|
726 _next->set(n, _first[*_highest]); |
|
727 if (_first[*_highest] != INVALID) { |
|
728 _prev->set(_first[*_highest], n); |
|
729 } else { |
|
730 _last[*_highest] = n; |
|
731 } |
|
732 _first[*_highest] = n; |
|
733 } |
|
734 } else { |
|
735 |
|
736 deactivate(n); |
|
737 if (!(*_active)[_first[*_highest]]) { |
|
738 ++_highest; |
|
739 if (_highest != _sets.back().end() && |
|
740 !(*_active)[_first[*_highest]]) { |
|
741 _highest = _sets.back().end(); |
|
742 } |
|
743 } |
|
744 } |
|
745 } |
|
746 |
|
747 if ((*_excess)[target] < _min_cut) { |
|
748 _min_cut = (*_excess)[target]; |
|
749 for (NodeIt i(_graph); i != INVALID; ++i) { |
|
750 _min_cut_map->set(i, false); |
|
751 } |
|
752 for (std::list<int>::iterator it = _sets.back().begin(); |
|
753 it != _sets.back().end(); ++it) { |
|
754 Node n = _first[*it]; |
|
755 while (n != INVALID) { |
|
756 _min_cut_map->set(n, true); |
|
757 n = (*_next)[n]; |
|
758 } |
|
759 } |
|
760 } |
|
761 |
|
762 { |
|
763 Node new_target; |
|
764 if ((*_prev)[target] != INVALID) { |
|
765 _last[(*_bucket)[target]] = (*_prev)[target]; |
|
766 _next->set((*_prev)[target], INVALID); |
|
767 new_target = (*_prev)[target]; |
|
768 } else { |
|
769 _sets.back().pop_back(); |
|
770 if (_sets.back().empty()) { |
|
771 _sets.pop_back(); |
|
772 if (_sets.empty()) |
|
773 break; |
|
774 for (std::list<int>::iterator it = _sets.back().begin(); |
|
775 it != _sets.back().end(); ++it) { |
|
776 _dormant[*it] = false; |
|
777 } |
|
778 } |
|
779 new_target = _last[_sets.back().back()]; |
|
780 } |
|
781 |
|
782 _bucket->set(target, 0); |
|
783 |
|
784 _source_set->set(target, true); |
|
785 for (InEdgeIt e(_graph, target); e != INVALID; ++e) { |
|
786 Value rem = (*_capacity)[e] - (*_flow)[e]; |
|
787 if (!_tolerance.positive(rem)) continue; |
|
788 Node v = _graph.source(e); |
|
789 if (!(*_active)[v] && !(*_source_set)[v]) { |
|
790 activate(v); |
|
791 } |
|
792 _excess->set(v, (*_excess)[v] + rem); |
|
793 _flow->set(e, (*_capacity)[e]); |
|
794 } |
|
795 |
|
796 for (OutEdgeIt e(_graph, target); e != INVALID; ++e) { |
|
797 Value rem = (*_flow)[e]; |
|
798 if (!_tolerance.positive(rem)) continue; |
|
799 Node v = _graph.target(e); |
|
800 if (!(*_active)[v] && !(*_source_set)[v]) { |
|
801 activate(v); |
|
802 } |
|
803 _excess->set(v, (*_excess)[v] + rem); |
|
804 _flow->set(e, 0); |
|
805 } |
|
806 |
|
807 target = new_target; |
|
808 if ((*_active)[target]) { |
|
809 deactivate(target); |
|
810 } |
|
811 |
|
812 _highest = _sets.back().begin(); |
|
813 while (_highest != _sets.back().end() && |
|
814 !(*_active)[_first[*_highest]]) { |
|
815 ++_highest; |
|
816 } |
|
817 } |
|
818 } |
|
819 } |
419 |
820 |
420 public: |
821 public: |
421 |
822 |
422 /// \name Execution control |
823 /// \name Execution control |
423 /// The simplest way to execute the algorithm is to use |
824 /// The simplest way to execute the algorithm is to use |
433 /// |
834 /// |
434 /// Initializes the internal data structures. It creates |
835 /// Initializes the internal data structures. It creates |
435 /// the maps, residual graph adaptors and some bucket structures |
836 /// the maps, residual graph adaptors and some bucket structures |
436 /// for the algorithm. |
837 /// for the algorithm. |
437 void init() { |
838 void init() { |
438 init(NodeIt(*_graph)); |
839 init(NodeIt(_graph)); |
439 } |
840 } |
440 |
841 |
441 /// \brief Initializes the internal data structures. |
842 /// \brief Initializes the internal data structures. |
442 /// |
843 /// |
443 /// Initializes the internal data structures. It creates |
844 /// Initializes the internal data structures. It creates |
444 /// the maps, residual graph adaptor and some bucket structures |
845 /// the maps, residual graph adaptor and some bucket structures |
445 /// for the algorithm. Node \c source is used as the push-relabel |
846 /// for the algorithm. Node \c source is used as the push-relabel |
446 /// algorithm's source. |
847 /// algorithm's source. |
447 void init(const Node& source) { |
848 void init(const Node& source) { |
448 _source = source; |
849 _source = source; |
449 _node_num = countNodes(*_graph); |
850 |
|
851 _node_num = countNodes(_graph); |
|
852 |
|
853 _first.resize(_node_num); |
|
854 _last.resize(_node_num); |
450 |
855 |
451 _dormant.resize(_node_num); |
856 _dormant.resize(_node_num); |
452 _level_size.resize(_node_num, 0); |
857 |
453 _active_nodes.resize(_node_num); |
858 if (!_flow) { |
454 |
859 _flow = new FlowMap(_graph); |
455 if (!_preflow) { |
860 } |
456 _preflow = new FlowMap(*_graph); |
861 if (!_next) { |
457 } |
862 _next = new typename Graph::template NodeMap<Node>(_graph); |
458 if (!_wake) { |
863 } |
459 _wake = new WakeMap(*_graph); |
864 if (!_prev) { |
460 } |
865 _prev = new typename Graph::template NodeMap<Node>(_graph); |
461 if (!_dist) { |
866 } |
462 _dist = new DistMap(*_graph); |
867 if (!_active) { |
|
868 _active = new typename Graph::template NodeMap<bool>(_graph); |
|
869 } |
|
870 if (!_bucket) { |
|
871 _bucket = new typename Graph::template NodeMap<int>(_graph); |
463 } |
872 } |
464 if (!_excess) { |
873 if (!_excess) { |
465 _excess = new ExcessMap(*_graph); |
874 _excess = new ExcessMap(_graph); |
466 } |
875 } |
467 if (!_source_set) { |
876 if (!_source_set) { |
468 _source_set = new SourceSetMap(*_graph); |
877 _source_set = new SourceSetMap(_graph); |
469 } |
|
470 if (!_out_res_graph) { |
|
471 _out_res_graph = new OutResGraph(*_graph, *_capacity, *_preflow); |
|
472 } |
|
473 if (!_rev_graph) { |
|
474 _rev_graph = new RevGraph(*_graph); |
|
475 } |
|
476 if (!_in_res_graph) { |
|
477 _in_res_graph = new InResGraph(*_rev_graph, *_capacity, *_preflow); |
|
478 } |
878 } |
479 if (!_min_cut_map) { |
879 if (!_min_cut_map) { |
480 _min_cut_map = new MinCutMap(*_graph); |
880 _min_cut_map = new MinCutMap(_graph); |
481 } |
881 } |
482 |
882 |
483 _min_cut = std::numeric_limits<Value>::max(); |
883 _min_cut = std::numeric_limits<Value>::max(); |
484 } |
884 } |
485 |
885 |
486 |
886 |
487 /// \brief Calculates a minimum cut with \f$ source \f$ on the |
887 /// \brief Calculates a minimum cut with \f$ source \f$ on the |
488 /// source-side. |
888 /// source-side. |
489 /// |
889 /// |
|
890 /// Calculates a minimum cut with \f$ source \f$ on the |
|
891 /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source |
|
892 /// \in X \f$ and minimal out-degree). |
|
893 void calculateOut() { |
|
894 findMinCutOut(); |
|
895 } |
|
896 |
490 /// \brief Calculates a minimum cut with \f$ source \f$ on the |
897 /// \brief Calculates a minimum cut with \f$ source \f$ on the |
491 /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$ |
898 /// target-side. |
492 /// and minimal out-degree). |
|
493 void calculateOut() { |
|
494 for (NodeIt it(*_graph); it != INVALID; ++it) { |
|
495 if (it != _source) { |
|
496 calculateOut(it); |
|
497 return; |
|
498 } |
|
499 } |
|
500 } |
|
501 |
|
502 /// \brief Calculates a minimum cut with \f$ source \f$ on the |
|
503 /// source-side. |
|
504 /// |
899 /// |
505 /// \brief Calculates a minimum cut with \f$ source \f$ on the |
900 /// Calculates a minimum cut with \f$ source \f$ on the |
506 /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$ |
901 /// target-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source |
507 /// and minimal out-degree). The \c target is the initial target |
902 /// \in X \f$ and minimal out-degree). |
508 /// for the push-relabel algorithm. |
|
509 void calculateOut(const Node& target) { |
|
510 findMinCut(target, true, *_out_res_graph); |
|
511 } |
|
512 |
|
513 /// \brief Calculates a minimum cut with \f$ source \f$ on the |
|
514 /// sink-side. |
|
515 /// |
|
516 /// \brief Calculates a minimum cut with \f$ source \f$ on the |
|
517 /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with |
|
518 /// \f$ source \notin X \f$ |
|
519 /// and minimal out-degree). |
|
520 void calculateIn() { |
903 void calculateIn() { |
521 for (NodeIt it(*_graph); it != INVALID; ++it) { |
904 findMinCutIn(); |
522 if (it != _source) { |
905 } |
523 calculateIn(it); |
906 |
524 return; |
|
525 } |
|
526 } |
|
527 } |
|
528 |
|
529 /// \brief Calculates a minimum cut with \f$ source \f$ on the |
|
530 /// sink-side. |
|
531 /// |
|
532 /// \brief Calculates a minimum cut with \f$ source \f$ on the |
|
533 /// sink-side (i.e. a set \f$ X\subsetneq V |
|
534 /// \f$ with \f$ source \notin X \f$ and minimal out-degree). |
|
535 /// The \c target is the initial |
|
536 /// target for the push-relabel algorithm. |
|
537 void calculateIn(const Node& target) { |
|
538 findMinCut(target, false, *_in_res_graph); |
|
539 } |
|
540 |
907 |
541 /// \brief Runs the algorithm. |
908 /// \brief Runs the algorithm. |
542 /// |
909 /// |
543 /// Runs the algorithm. It finds nodes \c source and \c target |
910 /// Runs the algorithm. It finds nodes \c source and \c target |
544 /// arbitrarily and then calls \ref init(), \ref calculateOut() |
911 /// arbitrarily and then calls \ref init(), \ref calculateOut() |
545 /// and \ref calculateIn(). |
912 /// and \ref calculateIn(). |
546 void run() { |
913 void run() { |
547 init(); |
914 init(); |
548 for (NodeIt it(*_graph); it != INVALID; ++it) { |
915 calculateOut(); |
549 if (it != _source) { |
916 calculateIn(); |
550 calculateOut(it); |
|
551 calculateIn(it); |
|
552 return; |
|
553 } |
|
554 } |
|
555 } |
917 } |
556 |
918 |
557 /// \brief Runs the algorithm. |
919 /// \brief Runs the algorithm. |
558 /// |
920 /// |
559 /// Runs the algorithm. It uses the given \c source node, finds a |
921 /// Runs the algorithm. It uses the given \c source node, finds a |
560 /// proper \c target and then calls the \ref init(), \ref |
922 /// proper \c target and then calls the \ref init(), \ref |
561 /// calculateOut() and \ref calculateIn(). |
923 /// calculateOut() and \ref calculateIn(). |
562 void run(const Node& s) { |
924 void run(const Node& s) { |
563 init(s); |
925 init(s); |
564 for (NodeIt it(*_graph); it != INVALID; ++it) { |
926 calculateOut(); |
565 if (it != _source) { |
927 calculateIn(); |
566 calculateOut(it); |
|
567 calculateIn(it); |
|
568 return; |
|
569 } |
|
570 } |
|
571 } |
|
572 |
|
573 /// \brief Runs the algorithm. |
|
574 /// |
|
575 /// Runs the algorithm. It just calls the \ref init() and then |
|
576 /// \ref calculateOut() and \ref calculateIn(). |
|
577 void run(const Node& s, const Node& t) { |
|
578 init(s); |
|
579 calculateOut(t); |
|
580 calculateIn(t); |
|
581 } |
928 } |
582 |
929 |
583 /// @} |
930 /// @} |
584 |
931 |
585 /// \name Query Functions |
932 /// \name Query Functions |