gtl/src/maxflow_pp.cpp

0001 /* This software is distributed under the GNU Lesser General Public License */
0002 //==========================================================================
0003 //
0004 //   maxflow_pp.cpp
0005 //
0006 //==========================================================================
0007 // $Id: maxflow_pp.cpp,v 1.7 2001/11/07 13:58:10 pick Exp $
0008
0009 #include <GTL/maxflow_pp.h>
0010
0011 #include <cstdlib>
0012 #include <iostream>
0013 #include <cassert>
0014
0015 #ifdef __GTL_MSVCC
0016 #   ifdef _DEBUG
0017 #   ifndef SEARCH_MEMORY_LEAKS_ENABLED
0018 #   error SEARCH NOT ENABLED
0019 #   endif
0020 #   define new DEBUG_NEW
0021 #   undef THIS_FILE
0022     static char THIS_FILE[] = __FILE__;
0023 #   endif   // _DEBUG
0024 #endif  // __GTL_MSVCC
0025
0026 __GTL_BEGIN_NAMESPACE
0027
0028 maxflow_pp::maxflow_pp()
0029 {
0030     max_graph_flow = 0.0;
0031     set_vars_executed = false;
0032 }
0033
0034
0035 maxflow_pp::~maxflow_pp()
0036 {
0037 }
0038
0039
0040 void maxflow_pp::set_vars(const edge_map<double>& edge_capacity)
0041 {
0042     this->edge_capacity = edge_capacity;
0043     artif_source_target = true;
0044     max_graph_flow = 0.0;
0045     set_vars_executed = true;
0046 }
0047
0048
0049 void maxflow_pp::set_vars(const edge_map<double>& edge_capacity,
0050     const node& net_source, const node& net_target)
0051 {
0052     this->edge_capacity = edge_capacity;
0053     this->net_source = net_source;
0054     this->net_target = net_target;
0055     artif_source_target = false;
0056     max_graph_flow = 0.0;
0057     set_vars_executed = true;
0058 }
0059
0060
0061 int maxflow_pp::check(graph& G)
0062 {
0063     if (!set_vars_executed)
0064     {
0065     return(GTL_ERROR);
0066     }
0067     graph::edge_iterator edge_it = G.edges_begin();
0068     graph::edge_iterator edges_end = G.edges_end();
0069     while (edge_it != edges_end)
0070     {
0071     if (edge_capacity[*edge_it] < 0)
0072     {
0073         return(GTL_ERROR);
0074     }
0075     ++edge_it;
0076     }
0077     // G.is_acyclic may be false
0078     if ((G.number_of_nodes() <= 1) || (!G.is_connected()) || (G.is_undirected()))
0079     {
0080     return(GTL_ERROR);
0081     }
0082     if (artif_source_target)
0083     {
0084     bool source_found = false;
0085     bool target_found = false;
0086     graph::node_iterator node_it = G.nodes_begin();
0087     graph::node_iterator nodes_end = G.nodes_end();
0088     while (node_it != nodes_end)
0089     {
0090         if ((*node_it).indeg() == 0)
0091         {
0092         source_found = true;
0093         }
0094         if ((*node_it).outdeg() == 0)
0095         {
0096         target_found = true;
0097         }
0098         ++node_it;
0099     }
0100     if (!(source_found && target_found))
0101     {
0102         return(GTL_ERROR);
0103     }
0104     }
0105     else
0106     {
0107     if (net_source == net_target)
0108     {
0109         return(GTL_ERROR);
0110     }
0111     }
0112     return(GTL_OK); // ok
0113 }
0114
0115
0116 int maxflow_pp::run(graph& G)
0117 {
0118     // init
0119     if (artif_source_target)
0120     {
0121     create_artif_source_target(G);
0122     }
0123     prepare_run(G);
0124
0125     double flow_value = 0;
0126     node min_tp_node;
0127     while(leveling(G) == TARGET_FROM_SOURCE_REACHABLE)
0128     {
0129     hide_unreachable_nodes(G);
0130     min_throughput_node(G, min_tp_node, flow_value);
0131     push(G, min_tp_node, flow_value);
0132     pull(G, min_tp_node, flow_value);
0133     comp_rem_net(G);
0134     }
0135
0136     restore_graph(G);
0137     return(GTL_OK);
0138 }
0139
0140
0141 double maxflow_pp::get_max_flow(const edge& e) const
0142 {
0143     return(edge_max_flow[e]);
0144 }
0145
0146
0147 double maxflow_pp::get_max_flow() const
0148 {
0149     return(max_graph_flow);
0150 }
0151
0152
0153 double maxflow_pp::get_rem_cap(const edge& e) const
0154 {
0155     return(edge_capacity[e] - edge_max_flow[e]);
0156 }
0157
0158
0159 void maxflow_pp::reset()
0160 {
0161 }
0162
0163
0164 void maxflow_pp::create_artif_source_target(graph& G)
0165 {
0166     net_source = G.new_node();
0167     net_target = G.new_node();
0168     edge e;
0169     graph::node_iterator node_it = G.nodes_begin();
0170     graph::node_iterator nodes_end = G.nodes_end();
0171     while (node_it != nodes_end)
0172     {
0173     if (*node_it != net_source && (*node_it).indeg() == 0)
0174     {
0175         e = G.new_edge(net_source, *node_it);
0176         edge_capacity[e] = 1.0; // 1.0 prevents e from hiding
0177         node::out_edges_iterator out_edge_it = (*node_it).out_edges_begin();
0178         node::out_edges_iterator out_edges_end = (*node_it).out_edges_end();
0179         while (out_edge_it != out_edges_end)
0180         {
0181         edge_capacity[e] += edge_capacity[*out_edge_it];
0182         ++out_edge_it;
0183         }
0184     }
0185     if (*node_it != net_target && (*node_it).outdeg() == 0)
0186     {
0187         e = G.new_edge(*node_it, net_target);
0188         edge_capacity[e] = 1.0; // 1.0 prevents e from hiding
0189         node::in_edges_iterator in_edge_it = (*node_it).in_edges_begin();
0190         node::in_edges_iterator in_edges_end = (*node_it).in_edges_end();
0191         while (in_edge_it != in_edges_end)
0192         {
0193         edge_capacity[e] += edge_capacity[*in_edge_it];
0194         ++in_edge_it;
0195         }
0196     }
0197     ++node_it;
0198     }
0199 }
0200
0201
0202 void maxflow_pp::prepare_run(const graph& G)
0203 {
0204     flow_update.init(G, 0.0);
0205     edge_max_flow.init(G, 0.0);
0206     edge_org.init(G, true);
0207     back_edge_exists.init(G, false);
0208     max_graph_flow = 0.0;
0209     full_edges.clear();
0210     temp_unvisible_nodes.clear();
0211     temp_unvisible_edges.clear();
0212 }
0213
0214
0215 int maxflow_pp::leveling(graph& G)
0216 {
0217     bool source_target_con = false;
0218     node_map<int> level(G, -1); // -1 means no level yet!
0219     queue<node> next_nodes;
0220     next_nodes.push(net_source);
0221     level[net_source] = 0;
0222     node cur_node;
0223     while (!next_nodes.empty())
0224     {
0225     cur_node = next_nodes.front();
0226     next_nodes.pop();
0227     node::out_edges_iterator out_edge_it = cur_node.out_edges_begin();
0228     node::out_edges_iterator out_edges_end = cur_node.out_edges_end();
0229     while (out_edge_it != out_edges_end)
0230     {
0231         if (level[(*out_edge_it).target()] == -1)
0232         {
0233         if ((*out_edge_it).target() == net_target)
0234         {
0235             source_target_con = true;
0236         }
0237         level[(*out_edge_it).target()] = level[cur_node] + 1;
0238         next_nodes.push((*out_edge_it).target());
0239         ++out_edge_it;
0240         }
0241         else if (level[(*out_edge_it).target()] <= level[cur_node])
0242         {
0243         node::out_edges_iterator temp_it = out_edge_it;
0244         ++out_edge_it;
0245         temp_unvisible_edges.push_back(*temp_it);
0246         G.hide_edge(*temp_it);
0247         }
0248         else
0249         {
0250         ++out_edge_it;
0251         }
0252     }
0253     }
0254     if (source_target_con)
0255     {
0256     return(TARGET_FROM_SOURCE_REACHABLE);
0257     }
0258     else
0259     {
0260     return(TARGET_FROM_SOURCE_NOT_REACHABLE);
0261     }
0262 }
0263
0264
0265 void maxflow_pp::hide_unreachable_nodes(graph& G)
0266 {
0267     node_map<bool> reachable_from_net_source(G, false);
0268     node_map<bool> reachable_from_net_target(G, false);
0269     queue<node> next_nodes;
0270     node cur_node;
0271
0272     next_nodes.push(net_source);
0273     reachable_from_net_source[net_source] = true;
0274     while (!next_nodes.empty())
0275     {
0276     cur_node = next_nodes.front();
0277     next_nodes.pop();
0278     node::out_edges_iterator out_edge_it = cur_node.out_edges_begin();
0279     node::out_edges_iterator out_edges_end = cur_node.out_edges_end();
0280     while (out_edge_it != out_edges_end)
0281     {
0282         node next = (*out_edge_it).target();
0283         if (!reachable_from_net_source[next])
0284         {
0285         next_nodes.push(next);
0286         reachable_from_net_source[next] = true;
0287         }
0288         ++out_edge_it;
0289     }
0290     }
0291
0292     next_nodes.push(net_target);
0293     reachable_from_net_target[net_target] = true;
0294     while (!next_nodes.empty())
0295     {
0296     cur_node = next_nodes.front();
0297     next_nodes.pop();
0298     node::in_edges_iterator in_edge_it = cur_node.in_edges_begin();
0299     node::in_edges_iterator in_edges_end = cur_node.in_edges_end();
0300     while (in_edge_it != in_edges_end)
0301     {
0302         node next = (*in_edge_it).source();
0303         if (!reachable_from_net_target[next])
0304         {
0305         next_nodes.push(next);
0306         reachable_from_net_target[next] = true;
0307         }
0308         ++in_edge_it;
0309     }
0310     }
0311
0312     graph::node_iterator node_it = G.nodes_begin();
0313     graph::node_iterator nodes_end = G.nodes_end();
0314     while (node_it != nodes_end)
0315     {
0316     if ((!reachable_from_net_source[*node_it]) ||
0317         (!reachable_from_net_target[*node_it]))
0318     {
0319         graph::node_iterator temp_it = node_it;
0320         ++node_it;
0321         temp_unvisible_nodes.push_back(*temp_it);
0322         store_temp_unvisible_edges(*temp_it);
0323         G.hide_node(*temp_it);
0324     }
0325     else
0326     {
0327         ++node_it;
0328     }
0329     }
0330 }
0331
0332
0333 void maxflow_pp::store_temp_unvisible_edges(const node& cur_node)
0334 {
0335     node::in_edges_iterator in_it = cur_node.in_edges_begin();
0336     node::in_edges_iterator in_edges_end = cur_node.in_edges_end();
0337     while (in_it != in_edges_end)
0338     {
0339     temp_unvisible_edges.push_back(*in_it);
0340     ++in_it;
0341     }
0342     node::out_edges_iterator out_it = cur_node.out_edges_begin();
0343     node::out_edges_iterator out_edges_end = cur_node.out_edges_end();
0344     while (out_it != out_edges_end)
0345     {
0346     temp_unvisible_edges.push_back(*out_it);
0347     ++out_it;
0348     }
0349 }
0350
0351
0352 void maxflow_pp::min_throughput_node(const graph& G, node& min_tp_node,
0353     double& flow_value)
0354 {
0355     min_tp_node = net_source;
0356     flow_value = comp_min_throughput(min_tp_node);
0357
0358     graph::node_iterator node_it = G.nodes_begin();
0359     graph::node_iterator nodes_end = G.nodes_end();
0360     double cur_tp;
0361     while (node_it != nodes_end)
0362     {
0363     cur_tp = comp_min_throughput(*node_it);
0364     if (cur_tp < flow_value)
0365     {
0366         min_tp_node = *node_it;
0367         flow_value = cur_tp;
0368     }
0369     ++node_it;
0370     }
0371 }
0372
0373
0374 double maxflow_pp::comp_min_throughput(const node cur_node) const
0375 {
0376     double in_flow = 0.0;
0377     double out_flow = 0.0;
0378     node::in_edges_iterator in_it = cur_node.in_edges_begin();
0379     node::in_edges_iterator in_edges_end = cur_node.in_edges_end();
0380     while (in_it != in_edges_end)
0381     {
0382     in_flow += edge_capacity[*in_it] - edge_max_flow[*in_it];
0383     ++in_it;
0384     }
0385     node::out_edges_iterator out_it = cur_node.out_edges_begin();
0386     node::out_edges_iterator out_edges_end = cur_node.out_edges_end();
0387     while (out_it != out_edges_end)
0388     {
0389     out_flow += edge_capacity[*out_it] - edge_max_flow[*out_it];
0390     ++out_it;
0391     }
0392     if (cur_node == net_source)
0393     {
0394     return(out_flow);
0395     }
0396     if (cur_node == net_target)
0397     {
0398     return(in_flow);
0399     }
0400     return(in_flow < out_flow ? in_flow : out_flow);
0401 }
0402
0403
0404 void maxflow_pp::get_sp_ahead(const graph& G, const node& start_node,
0405     node_map<edge>& last_edge)
0406 {
0407     queue<node> next_nodes;
0408     node_map<bool> visited(G, false);
0409     next_nodes.push(start_node);
0410     visited[start_node] = true;
0411
0412     node cur_node;
0413     while (!next_nodes.empty())
0414     {
0415     cur_node = next_nodes.front();
0416     next_nodes.pop();
0417
0418     node::out_edges_iterator out_edge_it = cur_node.out_edges_begin();
0419     node::out_edges_iterator out_edges_end = cur_node.out_edges_end();
0420     while (out_edge_it != out_edges_end)
0421     {
0422         node next = (*out_edge_it).target();
0423         if (!visited[next])
0424         {
0425         last_edge[next] = *out_edge_it;
0426         if (next == net_target)
0427         {
0428             return; // sp found
0429         }
0430         next_nodes.push(next);
0431         visited[next] = true;
0432         }
0433         ++out_edge_it;
0434     }
0435     }
0436 }
0437
0438
0439 void maxflow_pp::get_sp_backwards(const graph& G, const node& start_node,
0440     node_map<edge>& prev_edge)
0441 {
0442     queue<node> next_nodes;
0443     node_map<bool> visited(G, false);
0444     next_nodes.push(start_node);
0445     visited[start_node] = true;
0446
0447     node cur_node;
0448     while (!next_nodes.empty())
0449     {
0450     cur_node = next_nodes.front();
0451     next_nodes.pop();
0452
0453     node::in_edges_iterator in_edge_it = cur_node.in_edges_begin();
0454     node::in_edges_iterator in_edges_end = cur_node.in_edges_end();
0455     while (in_edge_it != in_edges_end)
0456     {
0457         node next = (*in_edge_it).source();
0458         if (!visited[next])
0459         {
0460         prev_edge[next] = *in_edge_it;
0461         if (next == net_source)
0462         {
0463             return; // sp found
0464         }
0465         next_nodes.push(next);
0466         visited[next] = true;
0467         }
0468         ++in_edge_it;
0469     }
0470     }
0471 }
0472
0473
0474 void maxflow_pp::push(graph& G, const node& start_node, const double flow_value)
0475 {
0476     node_map<edge> last_edge;
0477     double cur_flow = flow_value;
0478     double min_value = 0.0;
0479
0480     if (start_node == net_target)
0481     {
0482     return; // no push necessary
0483     }
0484     do
0485     {
0486     get_sp_ahead(G, start_node, last_edge);
0487     min_value = extra_charge_ahead(start_node, last_edge);
0488     if (min_value > cur_flow)
0489     {
0490         min_value = cur_flow;
0491     }
0492     node cur_node = net_target;
0493     do
0494     {
0495         if (edge_org[last_edge[cur_node]])
0496         {
0497         edge_max_flow[last_edge[cur_node]] += min_value;
0498         if (back_edge_exists[last_edge[cur_node]])
0499         {
0500             flow_update[back_edge[last_edge[cur_node]]] += min_value;
0501         }
0502         }
0503         else
0504         {
0505         edge_capacity[last_edge[cur_node]] -= min_value;
0506         flow_update[back_edge[last_edge[cur_node]]] += min_value;
0507         }
0508         if (edge_capacity[last_edge[cur_node]] <=
0509         edge_max_flow[last_edge[cur_node]])
0510         {
0511         full_edges.push_back(last_edge[cur_node]);
0512         G.hide_edge(last_edge[cur_node]);
0513         }
0514         cur_node = last_edge[cur_node].source();
0515     }
0516     while (cur_node != start_node);
0517     cur_flow -= min_value;
0518     if (cur_flow < 1e-015)  // quite hacky ;-)
0519     {
0520         cur_flow = 0.0; // to avoid rounding errors
0521     }
0522     } while (cur_flow > 0.0);
0523 }
0524
0525
0526 void maxflow_pp::pull(graph& G, const node& start_node, const double flow_value)
0527 {
0528     node_map<edge> prev_edge;
0529     double cur_flow = flow_value;
0530     double min_value = 0.0;
0531
0532     if (start_node == net_source)
0533     {
0534     return; // pull not necessary
0535     }
0536     do
0537     {
0538     get_sp_backwards(G, start_node, prev_edge);
0539     min_value = extra_charge_backwards(start_node, prev_edge);
0540     if (min_value > cur_flow)
0541     {
0542         min_value = cur_flow;
0543     }
0544     node cur_node = net_source;
0545     do
0546     {
0547         if (edge_org[prev_edge[cur_node]])
0548         {
0549         edge_max_flow[prev_edge[cur_node]] += min_value;
0550         if (back_edge_exists[prev_edge[cur_node]])
0551         {
0552             flow_update[back_edge[prev_edge[cur_node]]] += min_value;
0553         }
0554         }
0555         else
0556         {
0557         edge_capacity[prev_edge[cur_node]] -= min_value;
0558         flow_update[back_edge[prev_edge[cur_node]]] += min_value;
0559         }
0560         if (edge_capacity[prev_edge[cur_node]] <=
0561         edge_max_flow[prev_edge[cur_node]])
0562         {
0563         full_edges.push_back(prev_edge[cur_node]);
0564         G.hide_edge(prev_edge[cur_node]);
0565         }
0566         cur_node = prev_edge[cur_node].target();
0567     }
0568     while (cur_node != start_node);
0569     cur_flow -= min_value;
0570     if (cur_flow < 1e-015)  // quite hacky ;-)
0571     {
0572         cur_flow = 0.0; // to avoid rounding errors
0573     }
0574     } while (cur_flow > 0.0);
0575 }
0576
0577
0578 void maxflow_pp::comp_rem_net(graph& G)
0579 {
0580     // update back_edges
0581     graph::edge_iterator edge_it = G.edges_begin();
0582     graph::edge_iterator edges_end = G.edges_end();
0583     while (edge_it != edges_end)
0584     {
0585     single_edge_update(G, *edge_it);
0586     ++edge_it;
0587     }
0588     list<edge>::iterator list_it = full_edges.begin();
0589     list<edge>::iterator list_end = full_edges.end();
0590     while (list_it != list_end)
0591     {
0592     G.restore_edge(*list_it);
0593     if (flow_update[*list_it] > 0.0)
0594     {
0595         single_edge_update(G, *list_it);
0596         list<edge>::iterator temp_it = list_it;
0597         ++list_it;
0598         full_edges.erase(temp_it);  // now it's visible again
0599     }
0600     else
0601     {
0602         if (!back_edge_exists[*list_it])
0603         {
0604         create_back_edge(G, *list_it);
0605         edge_capacity[back_edge[*list_it]] = edge_max_flow[*list_it];
0606         }
0607         G.hide_edge(*list_it);
0608         ++list_it;
0609     }
0610     }
0611
0612
0613     // make hidden levels visible again
0614     list<node>::iterator temp_un_node_it = temp_unvisible_nodes.begin();
0615     list<node>::iterator temp_un_nodes_end = temp_unvisible_nodes.end();
0616     while (temp_un_node_it != temp_un_nodes_end)
0617     {
0618     G.restore_node(*temp_un_node_it);
0619     ++temp_un_node_it;
0620     }
0621     list<edge>::iterator temp_un_edge_it = temp_unvisible_edges.begin();
0622     list<edge>::iterator temp_un_edges_end = temp_unvisible_edges.end();
0623     while (temp_un_edge_it != temp_un_edges_end)
0624     {
0625     G.restore_edge(*temp_un_edge_it);
0626     if (flow_update[*temp_un_edge_it] > 0.0)
0627     {
0628         single_edge_update(G, *temp_un_edge_it);
0629     }
0630     ++temp_un_edge_it;
0631     }
0632     temp_unvisible_nodes.clear();
0633     temp_unvisible_edges.clear();
0634 }
0635
0636
0637 void maxflow_pp::single_edge_update(graph& G, edge cur_edge)
0638 {
0639     if(edge_org[cur_edge])
0640     {
0641     edge_max_flow[cur_edge] -= flow_update[cur_edge];
0642     flow_update[cur_edge] = 0.0;
0643     if (!back_edge_exists[cur_edge])
0644     {
0645         if (edge_max_flow[cur_edge] > 0.0)
0646         {
0647         create_back_edge(G, cur_edge);
0648         edge_capacity[back_edge[cur_edge]] = edge_max_flow[cur_edge];
0649         }
0650     }
0651     }
0652     else
0653     {
0654     edge_capacity[cur_edge] += flow_update[cur_edge];
0655     flow_update[cur_edge] = 0.0;
0656     }
0657 }
0658
0659
0660 double maxflow_pp::extra_charge_ahead(const node& start_node,
0661     const node_map<edge>& last_edge) const
0662 {
0663     node cur_node = net_target;
0664     double min_value = edge_capacity[last_edge[cur_node]]
0665     - edge_max_flow[last_edge[cur_node]];
0666     double cur_capacity;
0667
0668     do
0669     {
0670     cur_capacity = edge_capacity[last_edge[cur_node]]
0671         - edge_max_flow[last_edge[cur_node]];
0672     if (cur_capacity < min_value) min_value = cur_capacity;
0673     cur_node = last_edge[cur_node].source();
0674     }
0675     while (cur_node != start_node);
0676     return(min_value);
0677 }
0678
0679
0680 double maxflow_pp::extra_charge_backwards(const node& start_node, const node_map<edge>& prev_edge) const
0681 {
0682     node cur_node = net_source;
0683     double min_value = edge_capacity[prev_edge[cur_node]]
0684     - edge_max_flow[prev_edge[cur_node]];
0685     double cur_capacity;
0686
0687     do
0688     {
0689     cur_capacity = edge_capacity[prev_edge[cur_node]]
0690         - edge_max_flow[prev_edge[cur_node]];
0691     if (cur_capacity < min_value) min_value = cur_capacity;
0692     cur_node = prev_edge[cur_node].target();
0693     }
0694     while (cur_node != start_node);
0695     return(min_value);
0696 }
0697
0698
0699 void maxflow_pp::create_back_edge(graph& G, const edge& org_edge)
0700 {
0701     edge be = G.new_edge(org_edge.target(), org_edge.source());
0702     edge_org[be] = false;
0703     edges_not_org.push_back(be);
0704     back_edge[org_edge] = be;
0705     back_edge[be] = org_edge;
0706     edge_max_flow[be] = 0.0;
0707     edge_capacity[be] = 0.0;
0708     back_edge_exists[org_edge] = true;
0709     back_edge_exists[be] = true;    // a back edge always has a org. edge ;-)
0710     flow_update[be] = 0.0;
0711 }
0712
0713
0714 void maxflow_pp::comp_max_flow(const graph& G)
0715 {
0716     max_graph_flow = 0.0;
0717
0718     node::out_edges_iterator out_edge_it = net_source.out_edges_begin();
0719     node::out_edges_iterator out_edges_end = net_source.out_edges_end();
0720     while (out_edge_it != out_edges_end)
0721     {
0722     max_graph_flow += edge_max_flow[*out_edge_it];
0723     ++out_edge_it;
0724     }
0725 }
0726
0727
0728 void maxflow_pp::restore_graph(graph& G)
0729 {
0730     G.restore_graph();
0731     while (!edges_not_org.empty())
0732     {
0733     G.del_edge(edges_not_org.front());
0734     edges_not_org.pop_front();
0735     }
0736     comp_max_flow(G);
0737     if (artif_source_target)
0738     {
0739     G.del_node(net_source);
0740     G.del_node(net_target);
0741     }
0742 }
0743
0744 __GTL_END_NAMESPACE
0745
0746 //--------------------------------------------------------------------------
0747 //   end of file
0748 //--------------------------------------------------------------------------