// Copyright (C) 2001 Vladimir Prus <ghost@cs.msu.su>

// Copyright (C) 2001 Jeremy Siek <jsiek@cs.indiana.edu>

// 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)

// NOTE: this final is generated by libs/graph/doc/transitive_closure.w

#ifndef BOOST_GRAPH_TRANSITIVE_CLOSURE_HPP

#define BOOST_GRAPH_TRANSITIVE_CLOSURE_HPP

#include <vector>

#include <algorithm> // for std::min and std::max

#include <functional>

#include <boost/config.hpp>

#include <boost/bind.hpp>

#include <boost/graph/vector_as_graph.hpp>

#include <boost/graph/strong_components.hpp>

#include <boost/graph/topological_sort.hpp>

#include <boost/graph/graph_concepts.hpp>

#include <boost/graph/named_function_params.hpp>

namespace boost

{

  namespace detail

  {

    inline void

      union_successor_sets(const std::vector < std::size_t > &s1,

                           const std::vector < std::size_t > &s2,

                           std::vector < std::size_t > &s3)

    {

      BOOST_USING_STD_MIN();

      for (std::size_t k = 0; k < s1.size(); ++k)

        s3[k] = min BOOST_PREVENT_MACRO_SUBSTITUTION(s1[k], s2[k]);

    }

  }                             // namespace detail

  namespace detail

  {

    template < typename Container, typename ST = std::size_t,

      typename VT = typename Container::value_type >

      struct subscript_t:public std::unary_function < ST, VT >

    {

      typedef VT& result_type;

      subscript_t(Container & c):container(&c)

      {

      }

      VT & operator() (const ST & i) const

      {

        return (*container)[i];

      }

    protected:

        Container * container;

    };

    template < typename Container >

      subscript_t < Container > subscript(Container & c) {

      return subscript_t < Container > (c);

    }

  }                             // namespace detail

  template < typename Graph, typename GraphTC,

    typename G_to_TC_VertexMap,

    typename VertexIndexMap >

    void transitive_closure(const Graph & g, GraphTC & tc,

                            G_to_TC_VertexMap g_to_tc_map,

                            VertexIndexMap index_map)

  {

    if (num_vertices(g) == 0)

      return;

    typedef typename graph_traits < Graph >::vertex_descriptor vertex;

    typedef typename graph_traits < Graph >::edge_descriptor edge;

    typedef typename graph_traits < Graph >::vertex_iterator vertex_iterator;

    typedef typename property_traits < VertexIndexMap >::value_type size_type;

    typedef typename graph_traits <

      Graph >::adjacency_iterator adjacency_iterator;

    function_requires < VertexListGraphConcept < Graph > >();

    function_requires < AdjacencyGraphConcept < Graph > >();

    function_requires < VertexMutableGraphConcept < GraphTC > >();

    function_requires < EdgeMutableGraphConcept < GraphTC > >();

    function_requires < ReadablePropertyMapConcept < VertexIndexMap,

      vertex > >();

    typedef size_type cg_vertex;

    std::vector < cg_vertex > component_number_vec(num_vertices(g));

    iterator_property_map < cg_vertex *, VertexIndexMap, cg_vertex, cg_vertex& >

      component_number(&component_number_vec[0], index_map);

    int num_scc = strong_components(g, component_number,

                                    vertex_index_map(index_map));

    std::vector < std::vector < vertex > >components;

    build_component_lists(g, num_scc, component_number, components);

    typedef std::vector<std::vector<cg_vertex> > CG_t;

    CG_t CG(num_scc);

    for (cg_vertex s = 0; s < components.size(); ++s) {

      std::vector < cg_vertex > adj;

      for (size_type i = 0; i < components[s].size(); ++i) {

        vertex u = components[s][i];

        adjacency_iterator v, v_end;

        for (tie(v, v_end) = adjacent_vertices(u, g); v != v_end; ++v) {

          cg_vertex t = component_number[*v];

          if (s != t)           // Avoid loops in the condensation graph

            adj.push_back(t);

        }

      }

      std::sort(adj.begin(), adj.end());

      typename std::vector<cg_vertex>::iterator di =

        std::unique(adj.begin(), adj.end());

      if (di != adj.end())

        adj.erase(di, adj.end());

      CG[s] = adj;

    }

    std::vector<cg_vertex> topo_order;

    std::vector<cg_vertex> topo_number(num_vertices(CG));

    topological_sort(CG, std::back_inserter(topo_order),

                     vertex_index_map(identity_property_map()));

    std::reverse(topo_order.begin(), topo_order.end());

    size_type n = 0;

    for (typename std::vector<cg_vertex>::iterator iter = topo_order.begin();

         iter != topo_order.end(); ++iter)

      topo_number[*iter] = n++;

    for (size_type i = 0; i < num_vertices(CG); ++i)

      std::sort(CG[i].begin(), CG[i].end(),

                boost::bind(std::less<cg_vertex>(),

                            boost::bind(detail::subscript(topo_number), _1),

                            boost::bind(detail::subscript(topo_number), _2)));

    std::vector<std::vector<cg_vertex> > chains;

    {

      std::vector<cg_vertex> in_a_chain(num_vertices(CG));

      for (typename std::vector<cg_vertex>::iterator i = topo_order.begin();

           i != topo_order.end(); ++i) {

        cg_vertex v = *i;

        if (!in_a_chain[v]) {

          chains.resize(chains.size() + 1);

          std::vector<cg_vertex>& chain = chains.back();

          for (;;) {

            chain.push_back(v);

            in_a_chain[v] = true;

            typename graph_traits<CG_t>::adjacency_iterator adj_first, adj_last;

            tie(adj_first, adj_last) = adjacent_vertices(v, CG);

            typename graph_traits<CG_t>::adjacency_iterator next

              = std::find_if(adj_first, adj_last,

                             std::not1(detail::subscript(in_a_chain)));

            if (next != adj_last)

              v = *next;

            else

              break;            // end of chain, dead-end

          }

        }

      }

    }

    std::vector<size_type> chain_number(num_vertices(CG));

    std::vector<size_type> pos_in_chain(num_vertices(CG));

    for (size_type i = 0; i < chains.size(); ++i)

      for (size_type j = 0; j < chains[i].size(); ++j) {

        cg_vertex v = chains[i][j];

        chain_number[v] = i;

        pos_in_chain[v] = j;

      }

    cg_vertex inf = (std::numeric_limits< cg_vertex >::max)();

    std::vector<std::vector<cg_vertex> > successors(num_vertices(CG),

                                                    std::vector<cg_vertex>

                                                    (chains.size(), inf));

    for (typename std::vector<cg_vertex>::reverse_iterator

           i = topo_order.rbegin(); i != topo_order.rend(); ++i) {

      cg_vertex u = *i;

      typename graph_traits<CG_t>::adjacency_iterator adj, adj_last;

      for (tie(adj, adj_last) = adjacent_vertices(u, CG);

           adj != adj_last; ++adj) {

        cg_vertex v = *adj;

        if (topo_number[v] < successors[u][chain_number[v]]) {

          // Succ(u) = Succ(u) U Succ(v)

          detail::union_successor_sets(successors[u], successors[v],

                                       successors[u]);

          // Succ(u) = Succ(u) U {v}

          successors[u][chain_number[v]] = topo_number[v];

        }

      }

    }

    for (size_type i = 0; i < CG.size(); ++i)

      CG[i].clear();

    for (size_type i = 0; i < CG.size(); ++i)

      for (size_type j = 0; j < chains.size(); ++j) {

        size_type topo_num = successors[i][j];

        if (topo_num < inf) {

          cg_vertex v = topo_order[topo_num];

          for (size_type k = pos_in_chain[v]; k < chains[j].size(); ++k)

            CG[i].push_back(chains[j][k]);

        }

      }

    // Add vertices to the transitive closure graph

    typedef typename graph_traits < GraphTC >::vertex_descriptor tc_vertex;

    {

      vertex_iterator i, i_end;

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

        g_to_tc_map[*i] = add_vertex(tc);

    }

    // Add edges between all the vertices in two adjacent SCCs

    typename graph_traits<CG_t>::vertex_iterator si, si_end;

    for (tie(si, si_end) = vertices(CG); si != si_end; ++si) {

      cg_vertex s = *si;

      typename graph_traits<CG_t>::adjacency_iterator i, i_end;

      for (tie(i, i_end) = adjacent_vertices(s, CG); i != i_end; ++i) {

        cg_vertex t = *i;

        for (size_type k = 0; k < components[s].size(); ++k)

          for (size_type l = 0; l < components[t].size(); ++l)

            add_edge(g_to_tc_map[components[s][k]],

                     g_to_tc_map[components[t][l]], tc);

      }

    }

    // Add edges connecting all vertices in a SCC

    for (size_type i = 0; i < components.size(); ++i)

      if (components[i].size() > 1)

        for (size_type k = 0; k < components[i].size(); ++k)

          for (size_type l = 0; l < components[i].size(); ++l) {

            vertex u = components[i][k], v = components[i][l];

            add_edge(g_to_tc_map[u], g_to_tc_map[v], tc);

          }

    // Find loopbacks in the original graph.

    // Need to add it to transitive closure.

    {

      vertex_iterator i, i_end;

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

        {

          adjacency_iterator ab, ae;

          for (boost::tie(ab, ae) = adjacent_vertices(*i, g); ab != ae; ++ab)

            {

              if (*ab == *i)

                if (components[component_number[*i]].size() == 1)

                  add_edge(g_to_tc_map[*i], g_to_tc_map[*i], tc);

            }

        }

    }

  }

  template <typename Graph, typename GraphTC>

  void transitive_closure(const Graph & g, GraphTC & tc)

  {

    if (num_vertices(g) == 0)

      return;

    typedef typename property_map<Graph, vertex_index_t>::const_type

      VertexIndexMap;

    VertexIndexMap index_map = get(vertex_index, g);

    typedef typename graph_traits<GraphTC>::vertex_descriptor tc_vertex;

    std::vector<tc_vertex> to_tc_vec(num_vertices(g));

    iterator_property_map < tc_vertex *, VertexIndexMap, tc_vertex, tc_vertex&>

      g_to_tc_map(&to_tc_vec[0], index_map);

    transitive_closure(g, tc, g_to_tc_map, index_map);

  }

  namespace detail

  {

    template < typename Graph, typename GraphTC, typename G_to_TC_VertexMap,

      typename VertexIndexMap>

    void transitive_closure_dispatch

      (const Graph & g, GraphTC & tc,

       G_to_TC_VertexMap g_to_tc_map, VertexIndexMap index_map)

    {

      typedef typename graph_traits < GraphTC >::vertex_descriptor tc_vertex;

      typename std::vector < tc_vertex >::size_type

        n = is_default_param(g_to_tc_map) ? num_vertices(g) : 1;

      std::vector < tc_vertex > to_tc_vec(n);

      transitive_closure

        (g, tc,

         choose_param(g_to_tc_map, make_iterator_property_map

                      (to_tc_vec.begin(), index_map, to_tc_vec[0])),

         index_map);

    }

  }                             // namespace detail

  template < typename Graph, typename GraphTC,

    typename P, typename T, typename R >

    void transitive_closure(const Graph & g, GraphTC & tc,

                            const bgl_named_params < P, T, R > &params)

  {

    if (num_vertices(g) == 0)

      return;

    detail::transitive_closure_dispatch

      (g, tc, get_param(params, orig_to_copy_t()),

       choose_const_pmap(get_param(params, vertex_index), g, vertex_index) );

  }

  template < typename G > void warshall_transitive_closure(G & g)

  {

    typedef typename graph_traits < G >::vertex_descriptor vertex;

    typedef typename graph_traits < G >::vertex_iterator vertex_iterator;

    function_requires < AdjacencyMatrixConcept < G > >();

    function_requires < EdgeMutableGraphConcept < G > >();

    // Matrix form:

    // for k

    //  for i

    //    if A[i,k]

    //      for j

    //        A[i,j] = A[i,j] | A[k,j]

    vertex_iterator ki, ke, ii, ie, ji, je;

    for (tie(ki, ke) = vertices(g); ki != ke; ++ki)

      for (tie(ii, ie) = vertices(g); ii != ie; ++ii)

        if (edge(*ii, *ki, g).second)

          for (tie(ji, je) = vertices(g); ji != je; ++ji)

            if (!edge(*ii, *ji, g).second && edge(*ki, *ji, g).second) {

              add_edge(*ii, *ji, g);

            }

  }

  template < typename G > void warren_transitive_closure(G & g)

  {

    using namespace boost;

    typedef typename graph_traits < G >::vertex_descriptor vertex;

    typedef typename graph_traits < G >::vertex_iterator vertex_iterator;

    function_requires < AdjacencyMatrixConcept < G > >();

    function_requires < EdgeMutableGraphConcept < G > >();

    // Make sure second loop will work

    if (num_vertices(g) == 0)

      return;

    // for i = 2 to n

    //    for k = 1 to i - 1

    //      if A[i,k]

    //        for j = 1 to n

    //          A[i,j] = A[i,j] | A[k,j]

    vertex_iterator ic, ie, jc, je, kc, ke;

    for (tie(ic, ie) = vertices(g), ++ic; ic != ie; ++ic)

      for (tie(kc, ke) = vertices(g); *kc != *ic; ++kc)

        if (edge(*ic, *kc, g).second)

          for (tie(jc, je) = vertices(g); jc != je; ++jc)

            if (!edge(*ic, *jc, g).second && edge(*kc, *jc, g).second) {

              add_edge(*ic, *jc, g);

            }

    //  for i = 1 to n - 1

    //    for k = i + 1 to n

    //      if A[i,k]

    //        for j = 1 to n

    //          A[i,j] = A[i,j] | A[k,j]

    for (tie(ic, ie) = vertices(g), --ie; ic != ie; ++ic)

      for (kc = ic, ke = ie, ++kc; kc != ke; ++kc)

        if (edge(*ic, *kc, g).second)

          for (tie(jc, je) = vertices(g); jc != je; ++jc)

            if (!edge(*ic, *jc, g).second && edge(*kc, *jc, g).second) {

              add_edge(*ic, *jc, g);

            }

  }

}                               // namespace boost

#endif // BOOST_GRAPH_TRANSITIVE_CLOSURE_HPP