gtl/src/maxflow_sap.cpp

0001 /* This software is distributed under the GNU Lesser General Public License */
0002 //==========================================================================
0003 //
0004 //   maxflow_sap.cpp
0005 //
0006 //==========================================================================
0007 // $Id: maxflow_sap.cpp,v 1.6 2001/11/07 13:58:10 pick Exp $
0008
0009 #include <GTL/maxflow_sap.h>
0010 #include <cstdlib>
0011 #include <iostream>
0012 #include <cassert>
0013
0014 #ifdef __GTL_MSVCC
0015 #   ifdef _DEBUG
0016 #   ifndef SEARCH_MEMORY_LEAKS_ENABLED
0017 #   error SEARCH NOT ENABLED
0018 #   endif
0019 #   define new DEBUG_NEW
0020 #   undef THIS_FILE
0021     static char THIS_FILE[] = __FILE__;
0022 #   endif   // _DEBUG
0023 #endif  // __GTL_MSVCC
0024
0025 __GTL_BEGIN_NAMESPACE
0026
0027 maxflow_sap::maxflow_sap()
0028 {
0029     max_graph_flow = 0.0;
0030     set_vars_executed = false;
0031 }
0032
0033
0034 maxflow_sap::~maxflow_sap()
0035 {
0036 }
0037
0038
0039 void maxflow_sap::set_vars(const edge_map<double>& edge_capacity)
0040 {
0041     this->edge_capacity = edge_capacity;
0042     artif_source_target = true;
0043     max_graph_flow = 0.0;
0044     set_vars_executed = true;
0045 }
0046
0047
0048 void maxflow_sap::set_vars(const edge_map<double>& edge_capacity,
0049     const node& net_source, const node& net_target)
0050 {
0051     this->edge_capacity = edge_capacity;
0052     this->net_source = net_source;
0053     this->net_target = net_target;
0054     artif_source_target = false;
0055     max_graph_flow = 0.0;
0056     set_vars_executed = true;
0057 }
0058
0059
0060 int maxflow_sap::check(graph& G)
0061 {
0062     if (!set_vars_executed)
0063     {
0064         return(GTL_ERROR);
0065     }
0066     graph::edge_iterator edge_it = G.edges_begin();
0067     graph::edge_iterator edges_end = G.edges_end();
0068     while (edge_it != edges_end)
0069     {
0070         if (edge_capacity[*edge_it] < 0)
0071         {
0072             return(GTL_ERROR);
0073         }
0074         ++edge_it;
0075     }
0076     // G.is_acyclic may be false
0077     if ((G.number_of_nodes() <= 1) || (!G.is_connected())
0078         || (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); // everything ok
0113 }
0114
0115
0116 int maxflow_sap::run(graph& G)
0117 {
0118     // init
0119     if (artif_source_target)
0120     {
0121         create_artif_source_target(G);
0122     }
0123     bool go_on = true;
0124     node_map<edge> last_edge(G);
0125     node cur_node;
0126     int number_of_nodes = G.number_of_nodes();
0127     vector<int> numb(number_of_nodes, 0);
0128     prepare_run(G);
0129
0130     comp_dist_labels(G, numb);
0131     cur_node = net_source;
0132
0133     while (go_on)
0134     {
0135         if (has_an_admissible_arc(cur_node))
0136         {
0137             advance(cur_node, last_edge);
0138             if (cur_node == net_target)
0139             {
0140                 augment(G, last_edge);
0141                 cur_node = net_source;
0142             }
0143         }
0144         else    // only inadmissible edges
0145         {
0146             go_on = retreat(number_of_nodes, cur_node, last_edge, numb);
0147         }
0148     }
0149
0150     restore_graph(G);
0151     return(GTL_OK);
0152 }
0153
0154
0155 double maxflow_sap::get_max_flow(const edge& e) const
0156 {
0157     return(edge_max_flow[e]);
0158 }
0159
0160
0161 double maxflow_sap::get_max_flow() const
0162 {
0163     return(max_graph_flow);
0164 }
0165
0166
0167 double maxflow_sap::get_rem_cap(const edge& e) const
0168 {
0169     return(edge_capacity[e] - edge_max_flow[e]);
0170 }
0171
0172
0173 void maxflow_sap::reset()
0174 {
0175     max_graph_flow = 0.0;
0176     set_vars_executed = false;
0177 }
0178
0179
0180 void maxflow_sap::create_artif_source_target(graph& G)
0181 {
0182     net_source = G.new_node();
0183     net_target = G.new_node();
0184     edge e;
0185     graph::node_iterator node_it = G.nodes_begin();
0186     graph::node_iterator nodes_end = G.nodes_end();
0187     while (node_it != nodes_end)
0188     {
0189         if (*node_it != net_source && (*node_it).indeg() == 0)
0190         {
0191             e = G.new_edge(net_source, *node_it);
0192             edge_capacity[e] = 1.0; // 1.0 prevents e from hiding
0193             node::out_edges_iterator out_edge_it =
0194                 (*node_it).out_edges_begin();
0195             node::out_edges_iterator out_edges_end =
0196                 (*node_it).out_edges_end();
0197             while (out_edge_it != out_edges_end)
0198             {
0199                 edge_capacity[e] += edge_capacity[*out_edge_it];
0200                 ++out_edge_it;
0201             }
0202         }
0203         if (*node_it != net_target && (*node_it).outdeg() == 0)
0204         {
0205             e = G.new_edge(*node_it, net_target);
0206             edge_capacity[e] = 1.0; // 1.0 prevents e from hiding
0207             node::in_edges_iterator in_edge_it =
0208                 (*node_it).in_edges_begin();
0209             node::in_edges_iterator in_edges_end =
0210                 (*node_it).in_edges_end();
0211             while (in_edge_it != in_edges_end)
0212             {
0213                 edge_capacity[e] += edge_capacity[*in_edge_it];
0214                 ++in_edge_it;
0215             }
0216         }
0217         ++node_it;
0218     }
0219 }
0220
0221
0222 void maxflow_sap::prepare_run(const graph& G)
0223 {
0224     edge_max_flow.init(G, 0.0);
0225     edge_org.init(G, true);
0226     back_edge_exists.init(G, false);
0227     max_graph_flow = 0.0;
0228 }
0229
0230
0231 void maxflow_sap::comp_dist_labels(const graph& G, vector<int>& numb)
0232 {
0233     queue<node> next_nodes;
0234     node_map<bool> visited(G, false);
0235     next_nodes.push(net_target);
0236     visited[net_target] = true;
0237     dist_label[net_target] = 0;
0238     numb[0] = 1;    // only one sink
0239     node cur_node;
0240
0241     while (!next_nodes.empty())
0242     {
0243         cur_node = next_nodes.front();
0244         next_nodes.pop();
0245         node::in_edges_iterator in_edge_it = cur_node.in_edges_begin();
0246         node::in_edges_iterator in_edges_end = cur_node.in_edges_end();
0247         while (in_edge_it != in_edges_end)
0248         {
0249             node next = (*in_edge_it).source();
0250             if (!visited[next])
0251             {
0252                 next_nodes.push(next);
0253                 visited[next] = true;
0254                 dist_label[next] = dist_label[cur_node] + 1;
0255                 ++numb[dist_label[next]];
0256             }
0257             ++in_edge_it;
0258         }
0259     }
0260 }
0261
0262
0263 bool maxflow_sap::has_an_admissible_arc(const node cur_node)
0264 {
0265     node::out_edges_iterator out_edge_it = cur_node.out_edges_begin();
0266     node::out_edges_iterator out_edges_end = cur_node.out_edges_end();
0267     while (out_edge_it != out_edges_end)
0268     {
0269         if (dist_label[cur_node] == dist_label[(*out_edge_it).target()] + 1)
0270         {
0271             return true;
0272         }
0273         ++out_edge_it;
0274     }
0275     return false;
0276 }
0277
0278
0279 void maxflow_sap::advance(node& cur_node, node_map<edge>& last_edge)
0280 {
0281     node::out_edges_iterator out_edge_it = cur_node.out_edges_begin();
0282     node::out_edges_iterator out_edges_end = cur_node.out_edges_end();
0283     while (out_edge_it != out_edges_end)
0284     {
0285         if (dist_label[cur_node] == dist_label[(*out_edge_it).target()] + 1)
0286         {
0287             last_edge[(*out_edge_it).target()] = *out_edge_it;
0288             cur_node = (*out_edge_it).target();
0289         }
0290         ++out_edge_it;
0291     }
0292 }
0293
0294
0295 void maxflow_sap::augment(graph& G, const node_map<edge>& last_edge)
0296 {
0297     double additional_flow = free_capacity(last_edge);
0298     node cur_node = net_target;
0299
0300     do
0301     {
0302         if (edge_org[last_edge[cur_node]])
0303             // shortest path runs over a org. edge
0304         {
0305             if (!back_edge_exists[last_edge[cur_node]]) // create back edge
0306             {
0307                 create_back_edge(G, last_edge[cur_node]);
0308             }
0309             edge_max_flow[last_edge[cur_node]] += additional_flow;
0310             G.restore_edge(back_edge[last_edge[cur_node]]);
0311             edge_capacity[back_edge[last_edge[cur_node]]] +=
0312                 additional_flow;
0313         }
0314         else    // shortest path runs over a inserted back edge
0315         {
0316             edge oe = back_edge[last_edge[cur_node]];
0317             G.restore_edge(oe);
0318             edge_max_flow[oe] -= additional_flow;
0319             edge_capacity[last_edge[cur_node]] -= additional_flow;
0320         }
0321         if (edge_capacity[last_edge[cur_node]] <=
0322             edge_max_flow[last_edge[cur_node]])
0323         {
0324             G.hide_edge(last_edge[cur_node]);
0325         }
0326         cur_node = last_edge[cur_node].source();
0327     }
0328     while (cur_node != net_source);
0329 }
0330
0331
0332 bool maxflow_sap::retreat(const int number_of_nodes,
0333                           node& cur_node,
0334                           const node_map<edge>& last_edge,
0335                           vector<int>& numb)
0336 {
0337     --numb[dist_label[cur_node]];
0338     if (numb[dist_label[cur_node]] == 0)
0339     {
0340         return false;
0341     }
0342     else
0343     {
0344         dist_label[cur_node] =
0345             min_neighbour_label(number_of_nodes, cur_node) + 1;
0346         ++numb[dist_label[cur_node]];
0347         if (cur_node != net_source)
0348         {
0349             cur_node = last_edge[cur_node].source();
0350         }
0351         return true;
0352     }
0353 }
0354
0355
0356 int maxflow_sap::min_neighbour_label(const int number_of_nodes,
0357                                      const node cur_node) const
0358 {
0359     int min_value = number_of_nodes;    // if no out edge exists
0360
0361     node::out_edges_iterator out_edge_it = cur_node.out_edges_begin();
0362     node::out_edges_iterator out_edges_end = cur_node.out_edges_end();
0363     while (out_edge_it != out_edges_end)
0364     {
0365         if (min_value > dist_label[(*out_edge_it).target()])
0366         {
0367             min_value = dist_label[(*out_edge_it).target()];
0368         }
0369         ++out_edge_it;
0370     }
0371     return min_value;
0372 }
0373
0374
0375 double maxflow_sap::free_capacity(const node_map<edge>& last_edge) const
0376 {
0377     node cur_node = net_target;
0378     double min_value =
0379         edge_capacity[last_edge[cur_node]] -
0380         edge_max_flow[last_edge[cur_node]];
0381     double cur_capacity;
0382
0383     do
0384     {
0385         cur_capacity =
0386             edge_capacity[last_edge[cur_node]] -
0387             edge_max_flow[last_edge[cur_node]];
0388
0389         if (cur_capacity < min_value)
0390         {
0391             min_value = cur_capacity;
0392         }
0393         cur_node = last_edge[cur_node].source();
0394     }
0395     while (cur_node != net_source);
0396     return(min_value);
0397 }
0398
0399
0400 void maxflow_sap::create_back_edge(graph& G, const edge& org_edge)
0401 {
0402     edge be = G.new_edge(org_edge.target(), org_edge.source());
0403     edge_org[be] = false;
0404     edges_not_org.push_back(be);
0405     back_edge[org_edge] = be;
0406     back_edge[be] = org_edge;
0407     edge_max_flow[be] = 0.0;
0408     edge_capacity[be] = 0.0;
0409     back_edge_exists[org_edge] = true;
0410     back_edge_exists[be] = true;    // a back edge always has a org. edge
0411 }
0412
0413
0414 void maxflow_sap::comp_max_flow(const graph& G)
0415 {
0416     max_graph_flow = 0.0;
0417
0418     node::out_edges_iterator out_edge_it = net_source.out_edges_begin();
0419     node::out_edges_iterator out_edges_end = net_source.out_edges_end();
0420     while (out_edge_it != out_edges_end)
0421     {
0422         max_graph_flow += edge_max_flow[*out_edge_it];
0423         ++out_edge_it;
0424     }
0425 }
0426
0427
0428 void maxflow_sap::restore_graph(graph& G)
0429 {
0430     G.restore_graph();  // hidden edges can not be deleted!
0431     while (!edges_not_org.empty())
0432     {
0433         G.del_edge(edges_not_org.front());
0434         edges_not_org.pop_front();
0435     }
0436     comp_max_flow(G);
0437     if (artif_source_target)
0438     {
0439         G.del_node(net_source);
0440         G.del_node(net_target);
0441     }
0442 }
0443
0444 __GTL_END_NAMESPACE
0445
0446 //--------------------------------------------------------------------------
0447 //   end of file
0448 //--------------------------------------------------------------------------