// Copyright 2004 The Trustees of Indiana University.

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

 //  Authors: Douglas Gregor

 //           Andrew Lumsdaine

 #ifndef BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP

 #define BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP

 #include <cmath>

 #include <boost/graph/graph_traits.hpp>

 #include <boost/graph/named_function_params.hpp>

 #include <boost/graph/simple_point.hpp>

 #include <vector>

 #include <list>

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

 namespace boost {

 struct square_distance_attractive_force {

   template<typename Graph, typename T>

   T

   operator()(typename graph_traits<Graph>::edge_descriptor,

              T k,

              T d,

              const Graph&) const

   {

     return d * d / k;

   }

 };

 struct square_distance_repulsive_force {

   template<typename Graph, typename T>

   T

   operator()(typename graph_traits<Graph>::vertex_descriptor,

              typename graph_traits<Graph>::vertex_descriptor,

              T k,

              T d,

              const Graph&) const

   {

     return k * k / d;

   }

 };

 template<typename T>

 struct linear_cooling {

   typedef T result_type;

   linear_cooling(std::size_t iterations)

     : temp(T(iterations) / T(10)), step(0.1) { }

   linear_cooling(std::size_t iterations, T temp)

     : temp(temp), step(temp / T(iterations)) { }

   T operator()()

   {

     T old_temp = temp;

     temp -= step;

     if (temp < T(0)) temp = T(0);

     return old_temp;

   }

  private:

   T temp;

   T step;

 };

 struct all_force_pairs

 {

   template<typename Graph, typename ApplyForce >

   void operator()(const Graph& g, ApplyForce apply_force)

   {

     typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;

     vertex_iterator v, end;

     for (tie(v, end) = vertices(g); v != end; ++v) {

       vertex_iterator u = v;

       for (++u; u != end; ++u) {

         apply_force(*u, *v);

         apply_force(*v, *u);

       }

     }

   }

 };

 template<typename Dim, typename PositionMap>

 struct grid_force_pairs

 {

   template<typename Graph>

   explicit

   grid_force_pairs(Dim width, Dim height, PositionMap position, const Graph& g)

     : width(width), height(height), position(position)

   {

 #ifndef BOOST_NO_STDC_NAMESPACE

     using std::sqrt;

 #endif // BOOST_NO_STDC_NAMESPACE

     two_k = Dim(2) * sqrt(width*height / num_vertices(g));

   }

   template<typename Graph, typename ApplyForce >

   void operator()(const Graph& g, ApplyForce apply_force)

   {

     typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;

     typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;

     typedef std::list<vertex_descriptor> bucket_t;

     typedef std::vector<bucket_t> buckets_t;

     std::size_t columns = std::size_t(width / two_k + Dim(1));

     std::size_t rows = std::size_t(height / two_k + Dim(1));

     buckets_t buckets(rows * columns);

     vertex_iterator v, v_end;

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

       std::size_t column = std::size_t((position[*v].x + width  / 2) / two_k);

       std::size_t row    = std::size_t((position[*v].y + height / 2) / two_k);

       if (column >= columns) column = columns - 1;

       if (row >= rows) row = rows - 1;

       buckets[row * columns + column].push_back(*v);

     }

     for (std::size_t row = 0; row < rows; ++row)

       for (std::size_t column = 0; column < columns; ++column) {

         bucket_t& bucket = buckets[row * columns + column];

         typedef typename bucket_t::iterator bucket_iterator;

         for (bucket_iterator u = bucket.begin(); u != bucket.end(); ++u) {

           // Repulse vertices in this bucket

           bucket_iterator v = u;

           for (++v; v != bucket.end(); ++v) {

             apply_force(*u, *v);

             apply_force(*v, *u);

           }

           std::size_t adj_start_row = row == 0? 0 : row - 1;

           std::size_t adj_end_row = row == rows - 1? row : row + 1;

           std::size_t adj_start_column = column == 0? 0 : column - 1;

           std::size_t adj_end_column = column == columns - 1? column : column + 1;

           for (std::size_t other_row = adj_start_row; other_row <= adj_end_row;

                ++other_row)

             for (std::size_t other_column = adj_start_column; 

                  other_column <= adj_end_column; ++other_column)

               if (other_row != row || other_column != column) {

                 // Repulse vertices in this bucket

                 bucket_t& other_bucket 

                   = buckets[other_row * columns + other_column];

                 for (v = other_bucket.begin(); v != other_bucket.end(); ++v)

                   apply_force(*u, *v);

               }

         }

       }

   }

  private:

   Dim width;

   Dim height;

   PositionMap position;

   Dim two_k;

 };

 template<typename Dim, typename PositionMap, typename Graph>

 inline grid_force_pairs<Dim, PositionMap>

 make_grid_force_pairs(Dim width, Dim height, const PositionMap& position,

                       const Graph& g)

 { return grid_force_pairs<Dim, PositionMap>(width, height, position, g); }

 template<typename Graph, typename PositionMap, typename Dim>

 void

 scale_graph(const Graph& g, PositionMap position,

             Dim left, Dim top, Dim right, Dim bottom)

 {

   if (num_vertices(g) == 0) return;

   if (bottom > top) {

     using std::swap;

     swap(bottom, top);

   }

   typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;

   // Find min/max ranges

   Dim minX = position[*vertices(g).first].x, maxX = minX;

   Dim minY = position[*vertices(g).first].y, maxY = minY;

   vertex_iterator vi, vi_end;

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

     BOOST_USING_STD_MIN();

     BOOST_USING_STD_MAX();

     minX = min BOOST_PREVENT_MACRO_SUBSTITUTION (minX, position[*vi].x);

     maxX = max BOOST_PREVENT_MACRO_SUBSTITUTION (maxX, position[*vi].x);

     minY = min BOOST_PREVENT_MACRO_SUBSTITUTION (minY, position[*vi].y);

     maxY = max BOOST_PREVENT_MACRO_SUBSTITUTION (maxY, position[*vi].y);

   }

   // Scale to bounding box provided

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

     position[*vi].x = ((position[*vi].x - minX) / (maxX - minX))

                     * (right - left) + left;

     position[*vi].y = ((position[*vi].y - minY) / (maxY - minY))

                     * (top - bottom) + bottom;

   }

 }

 namespace detail {

   template<typename PositionMap, typename DisplacementMap,

            typename RepulsiveForce, typename Dim, typename Graph>

   struct fr_apply_force

   {

     typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;

     fr_apply_force(const PositionMap& position,

                    const DisplacementMap& displacement,

                    RepulsiveForce repulsive_force, Dim k, const Graph& g)

       : position(position), displacement(displacement),

         repulsive_force(repulsive_force), k(k), g(g)

     { }

     void operator()(vertex_descriptor u, vertex_descriptor v)

     {

 #ifndef BOOST_NO_STDC_NAMESPACE

       using std::sqrt;

 #endif // BOOST_NO_STDC_NAMESPACE

       if (u != v) {

         Dim delta_x = position[v].x - position[u].x;

         Dim delta_y = position[v].y - position[u].y;

         Dim dist = sqrt(delta_x * delta_x + delta_y * delta_y);

         Dim fr = repulsive_force(u, v, k, dist, g);

         displacement[v].x += delta_x / dist * fr;

         displacement[v].y += delta_y / dist * fr;

       }

     }

   private:

     PositionMap position;

     DisplacementMap displacement;

     RepulsiveForce repulsive_force;

     Dim k;

     const Graph& g;

   };

 } // end namespace detail

 template<typename Graph, typename PositionMap, typename Dim,

          typename AttractiveForce, typename RepulsiveForce,

          typename ForcePairs, typename Cooling, typename DisplacementMap>

 void

 fruchterman_reingold_force_directed_layout

  (const Graph&    g,

   PositionMap     position,

   Dim             width,

   Dim             height,

   AttractiveForce attractive_force,

   RepulsiveForce  repulsive_force,

   ForcePairs      force_pairs,

   Cooling         cool,

   DisplacementMap displacement)

 {

   typedef typename graph_traits<Graph>::vertex_iterator   vertex_iterator;

   typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;

   typedef typename graph_traits<Graph>::edge_iterator     edge_iterator;

 #ifndef BOOST_NO_STDC_NAMESPACE

   using std::sqrt;

 #endif // BOOST_NO_STDC_NAMESPACE

   Dim area = width * height;

   // assume positions are initialized randomly

   Dim k = sqrt(area / num_vertices(g));

   detail::fr_apply_force<PositionMap, DisplacementMap,

                          RepulsiveForce, Dim, Graph>

     apply_force(position, displacement, repulsive_force, k, g);

   Dim temp = cool();

   if (temp) do {

     // Calculate repulsive forces

     vertex_iterator v, v_end;

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

       displacement[*v].x = 0;

       displacement[*v].y = 0;

     }

     force_pairs(g, apply_force);

     // Calculate attractive forces

     edge_iterator e, e_end;

     for (tie(e, e_end) = edges(g); e != e_end; ++e) {

       vertex_descriptor v = source(*e, g);

       vertex_descriptor u = target(*e, g);

       Dim delta_x = position[v].x - position[u].x;

       Dim delta_y = position[v].y - position[u].y;

       Dim dist = sqrt(delta_x * delta_x + delta_y * delta_y);

       Dim fa = attractive_force(*e, k, dist, g);

       displacement[v].x -= delta_x / dist * fa;

       displacement[v].y -= delta_y / dist * fa;

       displacement[u].x += delta_x / dist * fa;

       displacement[u].y += delta_y / dist * fa;

     }

     // Update positions

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

       BOOST_USING_STD_MIN();

       BOOST_USING_STD_MAX();

       Dim disp_size = sqrt(displacement[*v].x * displacement[*v].x

                            + displacement[*v].y * displacement[*v].y);

       position[*v].x += displacement[*v].x / disp_size 

                      * min BOOST_PREVENT_MACRO_SUBSTITUTION (disp_size, temp);

       position[*v].y += displacement[*v].y / disp_size 

                      * min BOOST_PREVENT_MACRO_SUBSTITUTION (disp_size, temp);

       position[*v].x = min BOOST_PREVENT_MACRO_SUBSTITUTION 

                          (width / 2, 

                           max BOOST_PREVENT_MACRO_SUBSTITUTION(-width / 2, 

                                                                position[*v].x));

       position[*v].y = min BOOST_PREVENT_MACRO_SUBSTITUTION

                          (height / 2, 

                           max BOOST_PREVENT_MACRO_SUBSTITUTION(-height / 2, 

                                                                position[*v].y));

     }

   } while (temp = cool());

 }

 namespace detail {

   template<typename DisplacementMap>

   struct fr_force_directed_layout

   {

     template<typename Graph, typename PositionMap, typename Dim,

              typename AttractiveForce, typename RepulsiveForce,

              typename ForcePairs, typename Cooling,

              typename Param, typename Tag, typename Rest>

     static void

     run(const Graph&    g,

         PositionMap     position,

         Dim             width,

         Dim             height,

         AttractiveForce attractive_force,

         RepulsiveForce  repulsive_force,

         ForcePairs      force_pairs,

         Cooling         cool,

         DisplacementMap displacement,

         const bgl_named_params<Param, Tag, Rest>&)

     {

       fruchterman_reingold_force_directed_layout

         (g, position, width, height, attractive_force, repulsive_force,

          force_pairs, cool, displacement);

     }

   };

   template<>

   struct fr_force_directed_layout<error_property_not_found>

   {

     template<typename Graph, typename PositionMap, typename Dim,

              typename AttractiveForce, typename RepulsiveForce,

              typename ForcePairs, typename Cooling,

              typename Param, typename Tag, typename Rest>

     static void

     run(const Graph&    g,

         PositionMap     position,

         Dim             width,

         Dim             height,

         AttractiveForce attractive_force,

         RepulsiveForce  repulsive_force,

         ForcePairs      force_pairs,

         Cooling         cool,

         error_property_not_found,

         const bgl_named_params<Param, Tag, Rest>& params)

     {

       std::vector<simple_point<Dim> > displacements(num_vertices(g));

       fruchterman_reingold_force_directed_layout

         (g, position, width, height, attractive_force, repulsive_force,

          force_pairs, cool,

          make_iterator_property_map

          (displacements.begin(),

           choose_const_pmap(get_param(params, vertex_index), g,

                             vertex_index),

           simple_point<Dim>()));

     }

   };

 } // end namespace detail

 template<typename Graph, typename PositionMap, typename Dim, typename Param,

          typename Tag, typename Rest>

 void

 fruchterman_reingold_force_directed_layout

   (const Graph&    g,

    PositionMap     position,

    Dim             width,

    Dim             height,

    const bgl_named_params<Param, Tag, Rest>& params)

 {

   typedef typename property_value<bgl_named_params<Param,Tag,Rest>,

                                   vertex_displacement_t>::type D;

   detail::fr_force_directed_layout<D>::run

     (g, position, width, height,

      choose_param(get_param(params, attractive_force_t()),

                   square_distance_attractive_force()),

      choose_param(get_param(params, repulsive_force_t()),

                   square_distance_repulsive_force()),

      choose_param(get_param(params, force_pairs_t()),

                   make_grid_force_pairs(width, height, position, g)),

      choose_param(get_param(params, cooling_t()),

                   linear_cooling<Dim>(100)),

      get_param(params, vertex_displacement_t()),

      params);

 }

 template<typename Graph, typename PositionMap, typename Dim>

 void

 fruchterman_reingold_force_directed_layout(const Graph&    g,

                                            PositionMap     position,

                                            Dim             width,

                                            Dim             height)

 {

   fruchterman_reingold_force_directed_layout

     (g, position, width, height,

      attractive_force(square_distance_attractive_force()));

 }

 } // end namespace boost

 #endif // BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP