//=======================================================================

 // Copyright 2000 University of Notre Dame.

 // Authors: Jeremy G. Siek, Andrew Lumsdaine, Lie-Quan Lee

 //

 // Distributed under the Boost Software License, Version 1.0. (See

 // accompanying file LICENSE_1_0.txt or copy at

 // http://www.boost.org/LICENSE_1_0.txt)

 //=======================================================================

 #ifndef BOOST_EDGE_CONNECTIVITY

 #define BOOST_EDGE_CONNECTIVITY

 // WARNING: not-yet fully tested!

 #include <boost/config.hpp>

 #include <vector>

 #include <set>

 #include <algorithm>

 #include <boost/graph/edmunds_karp_max_flow.hpp>

 namespace boost {

   namespace detail {

     template <class Graph>

     inline

     std::pair<typename graph_traits<Graph>::vertex_descriptor,

               typename graph_traits<Graph>::degree_size_type>

     min_degree_vertex(Graph& g)

     {

       typedef graph_traits<Graph> Traits;

       typename Traits::vertex_descriptor p;

       typedef typename Traits::degree_size_type size_type;

       size_type delta = (std::numeric_limits<size_type>::max)();

       typename Traits::vertex_iterator i, iend;

       for (tie(i, iend) = vertices(g); i != iend; ++i)

         if (degree(*i, g) < delta) {

           delta = degree(*i, g);

           p = *i;

         }

       return std::make_pair(p, delta);

     }

     template <class Graph, class OutputIterator>

     void neighbors(const Graph& g, 

                    typename graph_traits<Graph>::vertex_descriptor u,

                    OutputIterator result)

     {

       typename graph_traits<Graph>::adjacency_iterator ai, aend;

       for (tie(ai, aend) = adjacent_vertices(u, g); ai != aend; ++ai)

         *result++ = *ai;

     }

     template <class Graph, class VertexIterator, class OutputIterator>

     void neighbors(const Graph& g, 

                    VertexIterator first, VertexIterator last,

                    OutputIterator result)

     {

       for (; first != last; ++first)

         neighbors(g, *first, result);

     }

   } // namespace detail

   // O(m n)

   template <class VertexListGraph, class OutputIterator>

   typename graph_traits<VertexListGraph>::degree_size_type

   edge_connectivity(VertexListGraph& g, OutputIterator disconnecting_set)

   {

     //-------------------------------------------------------------------------

     // Type Definitions

     typedef graph_traits<VertexListGraph> Traits;

     typedef typename Traits::vertex_iterator vertex_iterator;

     typedef typename Traits::edge_iterator edge_iterator;

     typedef typename Traits::out_edge_iterator out_edge_iterator;

     typedef typename Traits::vertex_descriptor vertex_descriptor;

     typedef typename Traits::degree_size_type degree_size_type;

     typedef color_traits<default_color_type> Color;

     typedef adjacency_list_traits<vecS, vecS, directedS> Tr;

     typedef typename Tr::edge_descriptor Tr_edge_desc;

     typedef adjacency_list<vecS, vecS, directedS, no_property, 

       property<edge_capacity_t, degree_size_type,

         property<edge_residual_capacity_t, degree_size_type,

           property<edge_reverse_t, Tr_edge_desc> > > > 

       FlowGraph;

     typedef typename graph_traits<FlowGraph>::edge_descriptor edge_descriptor;

     //-------------------------------------------------------------------------

     // Variable Declarations

     vertex_descriptor u, v, p, k;

     edge_descriptor e1, e2;

     bool inserted;

     vertex_iterator vi, vi_end;

     edge_iterator ei, ei_end;

     degree_size_type delta, alpha_star, alpha_S_k;

     std::set<vertex_descriptor> S, neighbor_S;

     std::vector<vertex_descriptor> S_star, non_neighbor_S;

     std::vector<default_color_type> color(num_vertices(g));

     std::vector<edge_descriptor> pred(num_vertices(g));

     //-------------------------------------------------------------------------

     // Create a network flow graph out of the undirected graph

     FlowGraph flow_g(num_vertices(g));

     typename property_map<FlowGraph, edge_capacity_t>::type

       cap = get(edge_capacity, flow_g);

     typename property_map<FlowGraph, edge_residual_capacity_t>::type

       res_cap = get(edge_residual_capacity, flow_g);

     typename property_map<FlowGraph, edge_reverse_t>::type

       rev_edge = get(edge_reverse, flow_g);

     for (tie(ei, ei_end) = edges(g); ei != ei_end; ++ei) {

       u = source(*ei, g), v = target(*ei, g);

       tie(e1, inserted) = add_edge(u, v, flow_g);

       cap[e1] = 1;

       tie(e2, inserted) = add_edge(v, u, flow_g);

       cap[e2] = 1; // not sure about this

       rev_edge[e1] = e2;

       rev_edge[e2] = e1;

     }

     //-------------------------------------------------------------------------

     // The Algorithm

     tie(p, delta) = detail::min_degree_vertex(g);

     S_star.push_back(p);

     alpha_star = delta;

     S.insert(p);

     neighbor_S.insert(p);

     detail::neighbors(g, S.begin(), S.end(), 

                       std::inserter(neighbor_S, neighbor_S.begin()));

     std::set_difference(vertices(g).first, vertices(g).second,

                         neighbor_S.begin(), neighbor_S.end(),

                         std::back_inserter(non_neighbor_S));

     while (!non_neighbor_S.empty()) { // at most n - 1 times

       k = non_neighbor_S.front();

       alpha_S_k = edmunds_karp_max_flow

         (flow_g, p, k, cap, res_cap, rev_edge, &color[0], &pred[0]);

       if (alpha_S_k < alpha_star) {

         alpha_star = alpha_S_k;

         S_star.clear();

         for (tie(vi, vi_end) = vertices(flow_g); vi != vi_end; ++vi)

           if (color[*vi] != Color::white())

             S_star.push_back(*vi);

       }

       S.insert(k);

       neighbor_S.insert(k);

       detail::neighbors(g, k, std::inserter(neighbor_S, neighbor_S.begin()));

       non_neighbor_S.clear();

       std::set_difference(vertices(g).first, vertices(g).second,

                           neighbor_S.begin(), neighbor_S.end(),

                           std::back_inserter(non_neighbor_S));

     }

     //-------------------------------------------------------------------------

     // Compute edges of the cut [S*, ~S*]

     std::vector<bool> in_S_star(num_vertices(g), false);

     typename std::vector<vertex_descriptor>::iterator si;

     for (si = S_star.begin(); si != S_star.end(); ++si)

       in_S_star[*si] = true;

     degree_size_type c = 0;

     for (si = S_star.begin(); si != S_star.end(); ++si) {

       out_edge_iterator ei, ei_end;

       for (tie(ei, ei_end) = out_edges(*si, g); ei != ei_end; ++ei)

         if (!in_S_star[target(*ei, g)]) {

           *disconnecting_set++ = *ei;

           ++c;

         }

     }

     return c;

   }

 } // namespace boost

 #endif // BOOST_EDGE_CONNECTIVITY