Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
maxflow_pp.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file maxflow_pp.cpp
1 /* This software is distributed under the GNU Lesser General Public License */
2 //==========================================================================
3 //
4 // maxflow_pp.cpp
5 //
6 //==========================================================================
7 // $Id: maxflow_pp.cpp,v 1.7 2001/11/07 13:58:10 pick Exp $
8 
9 #include <GTL/maxflow_pp.h>
10 
11 #include <cstdlib>
12 #include <iostream>
13 #include <cassert>
14 
15 #ifdef __GTL_MSVCC
16 # ifdef _DEBUG
17 # ifndef SEARCH_MEMORY_LEAKS_ENABLED
18 # error SEARCH NOT ENABLED
19 # endif
20 # define new DEBUG_NEW
21 # undef THIS_FILE
22  static char THIS_FILE[] = __FILE__;
23 # endif // _DEBUG
24 #endif // __GTL_MSVCC
25 
27 
29 {
30  max_graph_flow = 0.0;
31  set_vars_executed = false;
32 }
33 
34 
36 {
37 }
38 
39 
40 void maxflow_pp::set_vars(const edge_map<double>& edge_capacity)
41 {
42  this->edge_capacity = edge_capacity;
43  artif_source_target = true;
44  max_graph_flow = 0.0;
45  set_vars_executed = true;
46 }
47 
48 
49 void maxflow_pp::set_vars(const edge_map<double>& edge_capacity,
50  const node& net_source, const node& net_target)
51 {
52  this->edge_capacity = edge_capacity;
53  this->net_source = net_source;
54  this->net_target = net_target;
55  artif_source_target = false;
56  max_graph_flow = 0.0;
57  set_vars_executed = true;
58 }
59 
60 
62 {
63  if (!set_vars_executed)
64  {
65  return(GTL_ERROR);
66  }
67  graph::edge_iterator edge_it = G.edges_begin();
68  graph::edge_iterator edges_end = G.edges_end();
69  while (edge_it != edges_end)
70  {
71  if (edge_capacity[*edge_it] < 0)
72  {
73  return(GTL_ERROR);
74  }
75  ++edge_it;
76  }
77  // G.is_acyclic may be false
78  if ((G.number_of_nodes() <= 1) || (!G.is_connected()) || (G.is_undirected()))
79  {
80  return(GTL_ERROR);
81  }
83  {
84  bool source_found = false;
85  bool target_found = false;
86  graph::node_iterator node_it = G.nodes_begin();
87  graph::node_iterator nodes_end = G.nodes_end();
88  while (node_it != nodes_end)
89  {
90  if ((*node_it).indeg() == 0)
91  {
92  source_found = true;
93  }
94  if ((*node_it).outdeg() == 0)
95  {
96  target_found = true;
97  }
98  ++node_it;
99  }
100  if (!(source_found && target_found))
101  {
102  return(GTL_ERROR);
103  }
104  }
105  else
106  {
107  if (net_source == net_target)
108  {
109  return(GTL_ERROR);
110  }
111  }
112  return(GTL_OK); // ok
113 }
114 
115 
117 {
118  // init
120  {
122  }
123  prepare_run(G);
124 
125  double flow_value = 0;
126  node min_tp_node;
128  {
130  min_throughput_node(G, min_tp_node, flow_value);
131  push(G, min_tp_node, flow_value);
132  pull(G, min_tp_node, flow_value);
133  comp_rem_net(G);
134  }
135 
136  restore_graph(G);
137  return(GTL_OK);
138 }
139 
140 
141 double maxflow_pp::get_max_flow(const edge& e) const
142 {
143  return(edge_max_flow[e]);
144 }
145 
146 
148 {
149  return(max_graph_flow);
150 }
151 
152 
153 double maxflow_pp::get_rem_cap(const edge& e) const
154 {
155  return(edge_capacity[e] - edge_max_flow[e]);
156 }
157 
158 
160 {
161 }
162 
163 
165 {
166  net_source = G.new_node();
167  net_target = G.new_node();
168  edge e;
169  graph::node_iterator node_it = G.nodes_begin();
170  graph::node_iterator nodes_end = G.nodes_end();
171  while (node_it != nodes_end)
172  {
173  if (*node_it != net_source && (*node_it).indeg() == 0)
174  {
175  e = G.new_edge(net_source, *node_it);
176  edge_capacity[e] = 1.0; // 1.0 prevents e from hiding
177  node::out_edges_iterator out_edge_it = (*node_it).out_edges_begin();
178  node::out_edges_iterator out_edges_end = (*node_it).out_edges_end();
179  while (out_edge_it != out_edges_end)
180  {
181  edge_capacity[e] += edge_capacity[*out_edge_it];
182  ++out_edge_it;
183  }
184  }
185  if (*node_it != net_target && (*node_it).outdeg() == 0)
186  {
187  e = G.new_edge(*node_it, net_target);
188  edge_capacity[e] = 1.0; // 1.0 prevents e from hiding
189  node::in_edges_iterator in_edge_it = (*node_it).in_edges_begin();
190  node::in_edges_iterator in_edges_end = (*node_it).in_edges_end();
191  while (in_edge_it != in_edges_end)
192  {
193  edge_capacity[e] += edge_capacity[*in_edge_it];
194  ++in_edge_it;
195  }
196  }
197  ++node_it;
198  }
199 }
200 
201 
203 {
204  flow_update.init(G, 0.0);
205  edge_max_flow.init(G, 0.0);
206  edge_org.init(G, true);
207  back_edge_exists.init(G, false);
208  max_graph_flow = 0.0;
209  full_edges.clear();
210  temp_unvisible_nodes.clear();
211  temp_unvisible_edges.clear();
212 }
213 
214 
216 {
217  bool source_target_con = false;
218  node_map<int> level(G, -1); // -1 means no level yet!
219  queue<node> next_nodes;
220  next_nodes.push(net_source);
221  level[net_source] = 0;
222  node cur_node;
223  while (!next_nodes.empty())
224  {
225  cur_node = next_nodes.front();
226  next_nodes.pop();
227  node::out_edges_iterator out_edge_it = cur_node.out_edges_begin();
228  node::out_edges_iterator out_edges_end = cur_node.out_edges_end();
229  while (out_edge_it != out_edges_end)
230  {
231  if (level[(*out_edge_it).target()] == -1)
232  {
233  if ((*out_edge_it).target() == net_target)
234  {
235  source_target_con = true;
236  }
237  level[(*out_edge_it).target()] = level[cur_node] + 1;
238  next_nodes.push((*out_edge_it).target());
239  ++out_edge_it;
240  }
241  else if (level[(*out_edge_it).target()] <= level[cur_node])
242  {
243  node::out_edges_iterator temp_it = out_edge_it;
244  ++out_edge_it;
245  temp_unvisible_edges.push_back(*temp_it);
246  G.hide_edge(*temp_it);
247  }
248  else
249  {
250  ++out_edge_it;
251  }
252  }
253  }
254  if (source_target_con)
255  {
257  }
258  else
259  {
261  }
262 }
263 
264 
266 {
267  node_map<bool> reachable_from_net_source(G, false);
268  node_map<bool> reachable_from_net_target(G, false);
269  queue<node> next_nodes;
270  node cur_node;
271 
272  next_nodes.push(net_source);
273  reachable_from_net_source[net_source] = true;
274  while (!next_nodes.empty())
275  {
276  cur_node = next_nodes.front();
277  next_nodes.pop();
278  node::out_edges_iterator out_edge_it = cur_node.out_edges_begin();
279  node::out_edges_iterator out_edges_end = cur_node.out_edges_end();
280  while (out_edge_it != out_edges_end)
281  {
282  node next = (*out_edge_it).target();
283  if (!reachable_from_net_source[next])
284  {
285  next_nodes.push(next);
286  reachable_from_net_source[next] = true;
287  }
288  ++out_edge_it;
289  }
290  }
291 
292  next_nodes.push(net_target);
293  reachable_from_net_target[net_target] = true;
294  while (!next_nodes.empty())
295  {
296  cur_node = next_nodes.front();
297  next_nodes.pop();
298  node::in_edges_iterator in_edge_it = cur_node.in_edges_begin();
299  node::in_edges_iterator in_edges_end = cur_node.in_edges_end();
300  while (in_edge_it != in_edges_end)
301  {
302  node next = (*in_edge_it).source();
303  if (!reachable_from_net_target[next])
304  {
305  next_nodes.push(next);
306  reachable_from_net_target[next] = true;
307  }
308  ++in_edge_it;
309  }
310  }
311 
312  graph::node_iterator node_it = G.nodes_begin();
313  graph::node_iterator nodes_end = G.nodes_end();
314  while (node_it != nodes_end)
315  {
316  if ((!reachable_from_net_source[*node_it]) ||
317  (!reachable_from_net_target[*node_it]))
318  {
319  graph::node_iterator temp_it = node_it;
320  ++node_it;
321  temp_unvisible_nodes.push_back(*temp_it);
322  store_temp_unvisible_edges(*temp_it);
323  G.hide_node(*temp_it);
324  }
325  else
326  {
327  ++node_it;
328  }
329  }
330 }
331 
332 
334 {
335  node::in_edges_iterator in_it = cur_node.in_edges_begin();
336  node::in_edges_iterator in_edges_end = cur_node.in_edges_end();
337  while (in_it != in_edges_end)
338  {
339  temp_unvisible_edges.push_back(*in_it);
340  ++in_it;
341  }
342  node::out_edges_iterator out_it = cur_node.out_edges_begin();
343  node::out_edges_iterator out_edges_end = cur_node.out_edges_end();
344  while (out_it != out_edges_end)
345  {
346  temp_unvisible_edges.push_back(*out_it);
347  ++out_it;
348  }
349 }
350 
351 
352 void maxflow_pp::min_throughput_node(const graph& G, node& min_tp_node,
353  double& flow_value)
354 {
355  min_tp_node = net_source;
356  flow_value = comp_min_throughput(min_tp_node);
357 
358  graph::node_iterator node_it = G.nodes_begin();
359  graph::node_iterator nodes_end = G.nodes_end();
360  double cur_tp;
361  while (node_it != nodes_end)
362  {
363  cur_tp = comp_min_throughput(*node_it);
364  if (cur_tp < flow_value)
365  {
366  min_tp_node = *node_it;
367  flow_value = cur_tp;
368  }
369  ++node_it;
370  }
371 }
372 
373 
374 double maxflow_pp::comp_min_throughput(const node cur_node) const
375 {
376  double in_flow = 0.0;
377  double out_flow = 0.0;
378  node::in_edges_iterator in_it = cur_node.in_edges_begin();
379  node::in_edges_iterator in_edges_end = cur_node.in_edges_end();
380  while (in_it != in_edges_end)
381  {
382  in_flow += edge_capacity[*in_it] - edge_max_flow[*in_it];
383  ++in_it;
384  }
385  node::out_edges_iterator out_it = cur_node.out_edges_begin();
386  node::out_edges_iterator out_edges_end = cur_node.out_edges_end();
387  while (out_it != out_edges_end)
388  {
389  out_flow += edge_capacity[*out_it] - edge_max_flow[*out_it];
390  ++out_it;
391  }
392  if (cur_node == net_source)
393  {
394  return(out_flow);
395  }
396  if (cur_node == net_target)
397  {
398  return(in_flow);
399  }
400  return(in_flow < out_flow ? in_flow : out_flow);
401 }
402 
403 
404 void maxflow_pp::get_sp_ahead(const graph& G, const node& start_node,
405  node_map<edge>& last_edge)
406 {
407  queue<node> next_nodes;
408  node_map<bool> visited(G, false);
409  next_nodes.push(start_node);
410  visited[start_node] = true;
411 
412  node cur_node;
413  while (!next_nodes.empty())
414  {
415  cur_node = next_nodes.front();
416  next_nodes.pop();
417 
418  node::out_edges_iterator out_edge_it = cur_node.out_edges_begin();
419  node::out_edges_iterator out_edges_end = cur_node.out_edges_end();
420  while (out_edge_it != out_edges_end)
421  {
422  node next = (*out_edge_it).target();
423  if (!visited[next])
424  {
425  last_edge[next] = *out_edge_it;
426  if (next == net_target)
427  {
428  return; // sp found
429  }
430  next_nodes.push(next);
431  visited[next] = true;
432  }
433  ++out_edge_it;
434  }
435  }
436 }
437 
438 
439 void maxflow_pp::get_sp_backwards(const graph& G, const node& start_node,
440  node_map<edge>& prev_edge)
441 {
442  queue<node> next_nodes;
443  node_map<bool> visited(G, false);
444  next_nodes.push(start_node);
445  visited[start_node] = true;
446 
447  node cur_node;
448  while (!next_nodes.empty())
449  {
450  cur_node = next_nodes.front();
451  next_nodes.pop();
452 
453  node::in_edges_iterator in_edge_it = cur_node.in_edges_begin();
454  node::in_edges_iterator in_edges_end = cur_node.in_edges_end();
455  while (in_edge_it != in_edges_end)
456  {
457  node next = (*in_edge_it).source();
458  if (!visited[next])
459  {
460  prev_edge[next] = *in_edge_it;
461  if (next == net_source)
462  {
463  return; // sp found
464  }
465  next_nodes.push(next);
466  visited[next] = true;
467  }
468  ++in_edge_it;
469  }
470  }
471 }
472 
473 
474 void maxflow_pp::push(graph& G, const node& start_node, const double flow_value)
475 {
476  node_map<edge> last_edge;
477  double cur_flow = flow_value;
478  double min_value = 0.0;
479 
480  if (start_node == net_target)
481  {
482  return; // no push necessary
483  }
484  do
485  {
486  get_sp_ahead(G, start_node, last_edge);
487  min_value = extra_charge_ahead(start_node, last_edge);
488  if (min_value > cur_flow)
489  {
490  min_value = cur_flow;
491  }
492  node cur_node = net_target;
493  do
494  {
495  if (edge_org[last_edge[cur_node]])
496  {
497  edge_max_flow[last_edge[cur_node]] += min_value;
498  if (back_edge_exists[last_edge[cur_node]])
499  {
500  flow_update[back_edge[last_edge[cur_node]]] += min_value;
501  }
502  }
503  else
504  {
505  edge_capacity[last_edge[cur_node]] -= min_value;
506  flow_update[back_edge[last_edge[cur_node]]] += min_value;
507  }
508  if (edge_capacity[last_edge[cur_node]] <=
509  edge_max_flow[last_edge[cur_node]])
510  {
511  full_edges.push_back(last_edge[cur_node]);
512  G.hide_edge(last_edge[cur_node]);
513  }
514  cur_node = last_edge[cur_node].source();
515  }
516  while (cur_node != start_node);
517  cur_flow -= min_value;
518  if (cur_flow < 1e-015) // quite hacky ;-)
519  {
520  cur_flow = 0.0; // to avoid rounding errors
521  }
522  } while (cur_flow > 0.0);
523 }
524 
525 
526 void maxflow_pp::pull(graph& G, const node& start_node, const double flow_value)
527 {
528  node_map<edge> prev_edge;
529  double cur_flow = flow_value;
530  double min_value = 0.0;
531 
532  if (start_node == net_source)
533  {
534  return; // pull not necessary
535  }
536  do
537  {
538  get_sp_backwards(G, start_node, prev_edge);
539  min_value = extra_charge_backwards(start_node, prev_edge);
540  if (min_value > cur_flow)
541  {
542  min_value = cur_flow;
543  }
544  node cur_node = net_source;
545  do
546  {
547  if (edge_org[prev_edge[cur_node]])
548  {
549  edge_max_flow[prev_edge[cur_node]] += min_value;
550  if (back_edge_exists[prev_edge[cur_node]])
551  {
552  flow_update[back_edge[prev_edge[cur_node]]] += min_value;
553  }
554  }
555  else
556  {
557  edge_capacity[prev_edge[cur_node]] -= min_value;
558  flow_update[back_edge[prev_edge[cur_node]]] += min_value;
559  }
560  if (edge_capacity[prev_edge[cur_node]] <=
561  edge_max_flow[prev_edge[cur_node]])
562  {
563  full_edges.push_back(prev_edge[cur_node]);
564  G.hide_edge(prev_edge[cur_node]);
565  }
566  cur_node = prev_edge[cur_node].target();
567  }
568  while (cur_node != start_node);
569  cur_flow -= min_value;
570  if (cur_flow < 1e-015) // quite hacky ;-)
571  {
572  cur_flow = 0.0; // to avoid rounding errors
573  }
574  } while (cur_flow > 0.0);
575 }
576 
577 
579 {
580  // update back_edges
581  graph::edge_iterator edge_it = G.edges_begin();
582  graph::edge_iterator edges_end = G.edges_end();
583  while (edge_it != edges_end)
584  {
585  single_edge_update(G, *edge_it);
586  ++edge_it;
587  }
588  list<edge>::iterator list_it = full_edges.begin();
589  list<edge>::iterator list_end = full_edges.end();
590  while (list_it != list_end)
591  {
592  G.restore_edge(*list_it);
593  if (flow_update[*list_it] > 0.0)
594  {
595  single_edge_update(G, *list_it);
596  list<edge>::iterator temp_it = list_it;
597  ++list_it;
598  full_edges.erase(temp_it); // now it's visible again
599  }
600  else
601  {
602  if (!back_edge_exists[*list_it])
603  {
604  create_back_edge(G, *list_it);
605  edge_capacity[back_edge[*list_it]] = edge_max_flow[*list_it];
606  }
607  G.hide_edge(*list_it);
608  ++list_it;
609  }
610  }
611 
612 
613  // make hidden levels visible again
614  list<node>::iterator temp_un_node_it = temp_unvisible_nodes.begin();
615  list<node>::iterator temp_un_nodes_end = temp_unvisible_nodes.end();
616  while (temp_un_node_it != temp_un_nodes_end)
617  {
618  G.restore_node(*temp_un_node_it);
619  ++temp_un_node_it;
620  }
621  list<edge>::iterator temp_un_edge_it = temp_unvisible_edges.begin();
622  list<edge>::iterator temp_un_edges_end = temp_unvisible_edges.end();
623  while (temp_un_edge_it != temp_un_edges_end)
624  {
625  G.restore_edge(*temp_un_edge_it);
626  if (flow_update[*temp_un_edge_it] > 0.0)
627  {
628  single_edge_update(G, *temp_un_edge_it);
629  }
630  ++temp_un_edge_it;
631  }
632  temp_unvisible_nodes.clear();
633  temp_unvisible_edges.clear();
634 }
635 
636 
638 {
639  if(edge_org[cur_edge])
640  {
641  edge_max_flow[cur_edge] -= flow_update[cur_edge];
642  flow_update[cur_edge] = 0.0;
643  if (!back_edge_exists[cur_edge])
644  {
645  if (edge_max_flow[cur_edge] > 0.0)
646  {
647  create_back_edge(G, cur_edge);
648  edge_capacity[back_edge[cur_edge]] = edge_max_flow[cur_edge];
649  }
650  }
651  }
652  else
653  {
654  edge_capacity[cur_edge] += flow_update[cur_edge];
655  flow_update[cur_edge] = 0.0;
656  }
657 }
658 
659 
660 double maxflow_pp::extra_charge_ahead(const node& start_node,
661  const node_map<edge>& last_edge) const
662 {
663  node cur_node = net_target;
664  double min_value = edge_capacity[last_edge[cur_node]]
665  - edge_max_flow[last_edge[cur_node]];
666  double cur_capacity;
667 
668  do
669  {
670  cur_capacity = edge_capacity[last_edge[cur_node]]
671  - edge_max_flow[last_edge[cur_node]];
672  if (cur_capacity < min_value) min_value = cur_capacity;
673  cur_node = last_edge[cur_node].source();
674  }
675  while (cur_node != start_node);
676  return(min_value);
677 }
678 
679 
680 double maxflow_pp::extra_charge_backwards(const node& start_node, const node_map<edge>& prev_edge) const
681 {
682  node cur_node = net_source;
683  double min_value = edge_capacity[prev_edge[cur_node]]
684  - edge_max_flow[prev_edge[cur_node]];
685  double cur_capacity;
686 
687  do
688  {
689  cur_capacity = edge_capacity[prev_edge[cur_node]]
690  - edge_max_flow[prev_edge[cur_node]];
691  if (cur_capacity < min_value) min_value = cur_capacity;
692  cur_node = prev_edge[cur_node].target();
693  }
694  while (cur_node != start_node);
695  return(min_value);
696 }
697 
698 
699 void maxflow_pp::create_back_edge(graph& G, const edge& org_edge)
700 {
701  edge be = G.new_edge(org_edge.target(), org_edge.source());
702  edge_org[be] = false;
703  edges_not_org.push_back(be);
704  back_edge[org_edge] = be;
705  back_edge[be] = org_edge;
706  edge_max_flow[be] = 0.0;
707  edge_capacity[be] = 0.0;
708  back_edge_exists[org_edge] = true;
709  back_edge_exists[be] = true; // a back edge always has a org. edge ;-)
710  flow_update[be] = 0.0;
711 }
712 
713 
715 {
716  max_graph_flow = 0.0;
717 
720  while (out_edge_it != out_edges_end)
721  {
722  max_graph_flow += edge_max_flow[*out_edge_it];
723  ++out_edge_it;
724  }
725 }
726 
727 
729 {
730  G.restore_graph();
731  while (!edges_not_org.empty())
732  {
733  G.del_edge(edges_not_org.front());
734  edges_not_org.pop_front();
735  }
736  comp_max_flow(G);
738  {
739  G.del_node(net_source);
740  G.del_node(net_target);
741  }
742 }
743 
745 
746 //--------------------------------------------------------------------------
747 // end of file
748 //--------------------------------------------------------------------------