williamr@2: // Copyright 2004 The Trustees of Indiana University.
williamr@2: 
williamr@2: // Distributed under the Boost Software License, Version 1.0.
williamr@2: // (See accompanying file LICENSE_1_0.txt or copy at
williamr@2: // http://www.boost.org/LICENSE_1_0.txt)
williamr@2: 
williamr@2: //  Authors: Douglas Gregor
williamr@2: //           Andrew Lumsdaine
williamr@2: #ifndef BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP
williamr@2: #define BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP
williamr@2: 
williamr@2: #include <cmath>
williamr@2: #include <boost/graph/graph_traits.hpp>
williamr@2: #include <boost/graph/named_function_params.hpp>
williamr@2: #include <boost/graph/simple_point.hpp>
williamr@2: #include <vector>
williamr@2: #include <list>
williamr@2: #include <algorithm> // for std::min and std::max
williamr@2: 
williamr@2: namespace boost {
williamr@2: 
williamr@2: struct square_distance_attractive_force {
williamr@2:   template<typename Graph, typename T>
williamr@2:   T
williamr@2:   operator()(typename graph_traits<Graph>::edge_descriptor,
williamr@2:              T k,
williamr@2:              T d,
williamr@2:              const Graph&) const
williamr@2:   {
williamr@2:     return d * d / k;
williamr@2:   }
williamr@2: };
williamr@2: 
williamr@2: struct square_distance_repulsive_force {
williamr@2:   template<typename Graph, typename T>
williamr@2:   T
williamr@2:   operator()(typename graph_traits<Graph>::vertex_descriptor,
williamr@2:              typename graph_traits<Graph>::vertex_descriptor,
williamr@2:              T k,
williamr@2:              T d,
williamr@2:              const Graph&) const
williamr@2:   {
williamr@2:     return k * k / d;
williamr@2:   }
williamr@2: };
williamr@2: 
williamr@2: template<typename T>
williamr@2: struct linear_cooling {
williamr@2:   typedef T result_type;
williamr@2: 
williamr@2:   linear_cooling(std::size_t iterations)
williamr@2:     : temp(T(iterations) / T(10)), step(0.1) { }
williamr@2: 
williamr@2:   linear_cooling(std::size_t iterations, T temp)
williamr@2:     : temp(temp), step(temp / T(iterations)) { }
williamr@2: 
williamr@2:   T operator()()
williamr@2:   {
williamr@2:     T old_temp = temp;
williamr@2:     temp -= step;
williamr@2:     if (temp < T(0)) temp = T(0);
williamr@2:     return old_temp;
williamr@2:   }
williamr@2: 
williamr@2:  private:
williamr@2:   T temp;
williamr@2:   T step;
williamr@2: };
williamr@2: 
williamr@2: struct all_force_pairs
williamr@2: {
williamr@2:   template<typename Graph, typename ApplyForce >
williamr@2:   void operator()(const Graph& g, ApplyForce apply_force)
williamr@2:   {
williamr@2:     typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
williamr@2:     vertex_iterator v, end;
williamr@2:     for (tie(v, end) = vertices(g); v != end; ++v) {
williamr@2:       vertex_iterator u = v;
williamr@2:       for (++u; u != end; ++u) {
williamr@2:         apply_force(*u, *v);
williamr@2:         apply_force(*v, *u);
williamr@2:       }
williamr@2:     }
williamr@2:   }
williamr@2: };
williamr@2: 
williamr@2: template<typename Dim, typename PositionMap>
williamr@2: struct grid_force_pairs
williamr@2: {
williamr@2:   template<typename Graph>
williamr@2:   explicit
williamr@2:   grid_force_pairs(Dim width, Dim height, PositionMap position, const Graph& g)
williamr@2:     : width(width), height(height), position(position)
williamr@2:   {
williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE
williamr@2:     using std::sqrt;
williamr@2: #endif // BOOST_NO_STDC_NAMESPACE
williamr@2:     two_k = Dim(2) * sqrt(width*height / num_vertices(g));
williamr@2:   }
williamr@2: 
williamr@2:   template<typename Graph, typename ApplyForce >
williamr@2:   void operator()(const Graph& g, ApplyForce apply_force)
williamr@2:   {
williamr@2:     typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
williamr@2:     typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
williamr@2:     typedef std::list<vertex_descriptor> bucket_t;
williamr@2:     typedef std::vector<bucket_t> buckets_t;
williamr@2: 
williamr@2:     std::size_t columns = std::size_t(width / two_k + Dim(1));
williamr@2:     std::size_t rows = std::size_t(height / two_k + Dim(1));
williamr@2:     buckets_t buckets(rows * columns);
williamr@2:     vertex_iterator v, v_end;
williamr@2:     for (tie(v, v_end) = vertices(g); v != v_end; ++v) {
williamr@2:       std::size_t column = std::size_t((position[*v].x + width  / 2) / two_k);
williamr@2:       std::size_t row    = std::size_t((position[*v].y + height / 2) / two_k);
williamr@2: 
williamr@2:       if (column >= columns) column = columns - 1;
williamr@2:       if (row >= rows) row = rows - 1;
williamr@2:       buckets[row * columns + column].push_back(*v);
williamr@2:     }
williamr@2: 
williamr@2:     for (std::size_t row = 0; row < rows; ++row)
williamr@2:       for (std::size_t column = 0; column < columns; ++column) {
williamr@2:         bucket_t& bucket = buckets[row * columns + column];
williamr@2:         typedef typename bucket_t::iterator bucket_iterator;
williamr@2:         for (bucket_iterator u = bucket.begin(); u != bucket.end(); ++u) {
williamr@2:           // Repulse vertices in this bucket
williamr@2:           bucket_iterator v = u;
williamr@2:           for (++v; v != bucket.end(); ++v) {
williamr@2:             apply_force(*u, *v);
williamr@2:             apply_force(*v, *u);
williamr@2:           }
williamr@2: 
williamr@2:           std::size_t adj_start_row = row == 0? 0 : row - 1;
williamr@2:           std::size_t adj_end_row = row == rows - 1? row : row + 1;
williamr@2:           std::size_t adj_start_column = column == 0? 0 : column - 1;
williamr@2:           std::size_t adj_end_column = column == columns - 1? column : column + 1;
williamr@2:           for (std::size_t other_row = adj_start_row; other_row <= adj_end_row;
williamr@2:                ++other_row)
williamr@2:             for (std::size_t other_column = adj_start_column; 
williamr@2:                  other_column <= adj_end_column; ++other_column)
williamr@2:               if (other_row != row || other_column != column) {
williamr@2:                 // Repulse vertices in this bucket
williamr@2:                 bucket_t& other_bucket 
williamr@2:                   = buckets[other_row * columns + other_column];
williamr@2:                 for (v = other_bucket.begin(); v != other_bucket.end(); ++v)
williamr@2:                   apply_force(*u, *v);
williamr@2:               }
williamr@2:         }
williamr@2:       }
williamr@2:   }
williamr@2: 
williamr@2:  private:
williamr@2:   Dim width;
williamr@2:   Dim height;
williamr@2:   PositionMap position;
williamr@2:   Dim two_k;
williamr@2: };
williamr@2: 
williamr@2: template<typename Dim, typename PositionMap, typename Graph>
williamr@2: inline grid_force_pairs<Dim, PositionMap>
williamr@2: make_grid_force_pairs(Dim width, Dim height, const PositionMap& position,
williamr@2:                       const Graph& g)
williamr@2: { return grid_force_pairs<Dim, PositionMap>(width, height, position, g); }
williamr@2: 
williamr@2: template<typename Graph, typename PositionMap, typename Dim>
williamr@2: void
williamr@2: scale_graph(const Graph& g, PositionMap position,
williamr@2:             Dim left, Dim top, Dim right, Dim bottom)
williamr@2: {
williamr@2:   if (num_vertices(g) == 0) return;
williamr@2: 
williamr@2:   if (bottom > top) {
williamr@2:     using std::swap;
williamr@2:     swap(bottom, top);
williamr@2:   }
williamr@2: 
williamr@2:   typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
williamr@2: 
williamr@2:   // Find min/max ranges
williamr@2:   Dim minX = position[*vertices(g).first].x, maxX = minX;
williamr@2:   Dim minY = position[*vertices(g).first].y, maxY = minY;
williamr@2:   vertex_iterator vi, vi_end;
williamr@2:   for (tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) {
williamr@2:     BOOST_USING_STD_MIN();
williamr@2:     BOOST_USING_STD_MAX();
williamr@2:     minX = min BOOST_PREVENT_MACRO_SUBSTITUTION (minX, position[*vi].x);
williamr@2:     maxX = max BOOST_PREVENT_MACRO_SUBSTITUTION (maxX, position[*vi].x);
williamr@2:     minY = min BOOST_PREVENT_MACRO_SUBSTITUTION (minY, position[*vi].y);
williamr@2:     maxY = max BOOST_PREVENT_MACRO_SUBSTITUTION (maxY, position[*vi].y);
williamr@2:   }
williamr@2: 
williamr@2:   // Scale to bounding box provided
williamr@2:   for (tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) {
williamr@2:     position[*vi].x = ((position[*vi].x - minX) / (maxX - minX))
williamr@2:                     * (right - left) + left;
williamr@2:     position[*vi].y = ((position[*vi].y - minY) / (maxY - minY))
williamr@2:                     * (top - bottom) + bottom;
williamr@2:   }
williamr@2: }
williamr@2: 
williamr@2: namespace detail {
williamr@2:   template<typename PositionMap, typename DisplacementMap,
williamr@2:            typename RepulsiveForce, typename Dim, typename Graph>
williamr@2:   struct fr_apply_force
williamr@2:   {
williamr@2:     typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
williamr@2: 
williamr@2:     fr_apply_force(const PositionMap& position,
williamr@2:                    const DisplacementMap& displacement,
williamr@2:                    RepulsiveForce repulsive_force, Dim k, const Graph& g)
williamr@2:       : position(position), displacement(displacement),
williamr@2:         repulsive_force(repulsive_force), k(k), g(g)
williamr@2:     { }
williamr@2: 
williamr@2:     void operator()(vertex_descriptor u, vertex_descriptor v)
williamr@2:     {
williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE
williamr@2:       using std::sqrt;
williamr@2: #endif // BOOST_NO_STDC_NAMESPACE
williamr@2:       if (u != v) {
williamr@2:         Dim delta_x = position[v].x - position[u].x;
williamr@2:         Dim delta_y = position[v].y - position[u].y;
williamr@2:         Dim dist = sqrt(delta_x * delta_x + delta_y * delta_y);
williamr@2:         Dim fr = repulsive_force(u, v, k, dist, g);
williamr@2:         displacement[v].x += delta_x / dist * fr;
williamr@2:         displacement[v].y += delta_y / dist * fr;
williamr@2:       }
williamr@2:     }
williamr@2: 
williamr@2:   private:
williamr@2:     PositionMap position;
williamr@2:     DisplacementMap displacement;
williamr@2:     RepulsiveForce repulsive_force;
williamr@2:     Dim k;
williamr@2:     const Graph& g;
williamr@2:   };
williamr@2: 
williamr@2: } // end namespace detail
williamr@2: 
williamr@2: template<typename Graph, typename PositionMap, typename Dim,
williamr@2:          typename AttractiveForce, typename RepulsiveForce,
williamr@2:          typename ForcePairs, typename Cooling, typename DisplacementMap>
williamr@2: void
williamr@2: fruchterman_reingold_force_directed_layout
williamr@2:  (const Graph&    g,
williamr@2:   PositionMap     position,
williamr@2:   Dim             width,
williamr@2:   Dim             height,
williamr@2:   AttractiveForce attractive_force,
williamr@2:   RepulsiveForce  repulsive_force,
williamr@2:   ForcePairs      force_pairs,
williamr@2:   Cooling         cool,
williamr@2:   DisplacementMap displacement)
williamr@2: {
williamr@2:   typedef typename graph_traits<Graph>::vertex_iterator   vertex_iterator;
williamr@2:   typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
williamr@2:   typedef typename graph_traits<Graph>::edge_iterator     edge_iterator;
williamr@2: 
williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE
williamr@2:   using std::sqrt;
williamr@2: #endif // BOOST_NO_STDC_NAMESPACE
williamr@2: 
williamr@2:   Dim area = width * height;
williamr@2:   // assume positions are initialized randomly
williamr@2:   Dim k = sqrt(area / num_vertices(g));
williamr@2: 
williamr@2:   detail::fr_apply_force<PositionMap, DisplacementMap,
williamr@2:                          RepulsiveForce, Dim, Graph>
williamr@2:     apply_force(position, displacement, repulsive_force, k, g);
williamr@2: 
williamr@2:   Dim temp = cool();
williamr@2:   if (temp) do {
williamr@2:     // Calculate repulsive forces
williamr@2:     vertex_iterator v, v_end;
williamr@2:     for (tie(v, v_end) = vertices(g); v != v_end; ++v) {
williamr@2:       displacement[*v].x = 0;
williamr@2:       displacement[*v].y = 0;
williamr@2:     }
williamr@2:     force_pairs(g, apply_force);
williamr@2: 
williamr@2:     // Calculate attractive forces
williamr@2:     edge_iterator e, e_end;
williamr@2:     for (tie(e, e_end) = edges(g); e != e_end; ++e) {
williamr@2:       vertex_descriptor v = source(*e, g);
williamr@2:       vertex_descriptor u = target(*e, g);
williamr@2:       Dim delta_x = position[v].x - position[u].x;
williamr@2:       Dim delta_y = position[v].y - position[u].y;
williamr@2:       Dim dist = sqrt(delta_x * delta_x + delta_y * delta_y);
williamr@2:       Dim fa = attractive_force(*e, k, dist, g);
williamr@2: 
williamr@2:       displacement[v].x -= delta_x / dist * fa;
williamr@2:       displacement[v].y -= delta_y / dist * fa;
williamr@2:       displacement[u].x += delta_x / dist * fa;
williamr@2:       displacement[u].y += delta_y / dist * fa;
williamr@2:     }
williamr@2: 
williamr@2:     // Update positions
williamr@2:     for (tie(v, v_end) = vertices(g); v != v_end; ++v) {
williamr@2:       BOOST_USING_STD_MIN();
williamr@2:       BOOST_USING_STD_MAX();
williamr@2:       Dim disp_size = sqrt(displacement[*v].x * displacement[*v].x
williamr@2:                            + displacement[*v].y * displacement[*v].y);
williamr@2:       position[*v].x += displacement[*v].x / disp_size 
williamr@2:                      * min BOOST_PREVENT_MACRO_SUBSTITUTION (disp_size, temp);
williamr@2:       position[*v].y += displacement[*v].y / disp_size 
williamr@2:                      * min BOOST_PREVENT_MACRO_SUBSTITUTION (disp_size, temp);
williamr@2:       position[*v].x = min BOOST_PREVENT_MACRO_SUBSTITUTION 
williamr@2:                          (width / 2, 
williamr@2:                           max BOOST_PREVENT_MACRO_SUBSTITUTION(-width / 2, 
williamr@2:                                                                position[*v].x));
williamr@2:       position[*v].y = min BOOST_PREVENT_MACRO_SUBSTITUTION
williamr@2:                          (height / 2, 
williamr@2:                           max BOOST_PREVENT_MACRO_SUBSTITUTION(-height / 2, 
williamr@2:                                                                position[*v].y));
williamr@2:     }
williamr@2:   } while (temp = cool());
williamr@2: }
williamr@2: 
williamr@2: namespace detail {
williamr@2:   template<typename DisplacementMap>
williamr@2:   struct fr_force_directed_layout
williamr@2:   {
williamr@2:     template<typename Graph, typename PositionMap, typename Dim,
williamr@2:              typename AttractiveForce, typename RepulsiveForce,
williamr@2:              typename ForcePairs, typename Cooling,
williamr@2:              typename Param, typename Tag, typename Rest>
williamr@2:     static void
williamr@2:     run(const Graph&    g,
williamr@2:         PositionMap     position,
williamr@2:         Dim             width,
williamr@2:         Dim             height,
williamr@2:         AttractiveForce attractive_force,
williamr@2:         RepulsiveForce  repulsive_force,
williamr@2:         ForcePairs      force_pairs,
williamr@2:         Cooling         cool,
williamr@2:         DisplacementMap displacement,
williamr@2:         const bgl_named_params<Param, Tag, Rest>&)
williamr@2:     {
williamr@2:       fruchterman_reingold_force_directed_layout
williamr@2:         (g, position, width, height, attractive_force, repulsive_force,
williamr@2:          force_pairs, cool, displacement);
williamr@2:     }
williamr@2:   };
williamr@2: 
williamr@2:   template<>
williamr@2:   struct fr_force_directed_layout<error_property_not_found>
williamr@2:   {
williamr@2:     template<typename Graph, typename PositionMap, typename Dim,
williamr@2:              typename AttractiveForce, typename RepulsiveForce,
williamr@2:              typename ForcePairs, typename Cooling,
williamr@2:              typename Param, typename Tag, typename Rest>
williamr@2:     static void
williamr@2:     run(const Graph&    g,
williamr@2:         PositionMap     position,
williamr@2:         Dim             width,
williamr@2:         Dim             height,
williamr@2:         AttractiveForce attractive_force,
williamr@2:         RepulsiveForce  repulsive_force,
williamr@2:         ForcePairs      force_pairs,
williamr@2:         Cooling         cool,
williamr@2:         error_property_not_found,
williamr@2:         const bgl_named_params<Param, Tag, Rest>& params)
williamr@2:     {
williamr@2:       std::vector<simple_point<Dim> > displacements(num_vertices(g));
williamr@2:       fruchterman_reingold_force_directed_layout
williamr@2:         (g, position, width, height, attractive_force, repulsive_force,
williamr@2:          force_pairs, cool,
williamr@2:          make_iterator_property_map
williamr@2:          (displacements.begin(),
williamr@2:           choose_const_pmap(get_param(params, vertex_index), g,
williamr@2:                             vertex_index),
williamr@2:           simple_point<Dim>()));
williamr@2:     }
williamr@2:   };
williamr@2: 
williamr@2: } // end namespace detail
williamr@2: 
williamr@2: template<typename Graph, typename PositionMap, typename Dim, typename Param,
williamr@2:          typename Tag, typename Rest>
williamr@2: void
williamr@2: fruchterman_reingold_force_directed_layout
williamr@2:   (const Graph&    g,
williamr@2:    PositionMap     position,
williamr@2:    Dim             width,
williamr@2:    Dim             height,
williamr@2:    const bgl_named_params<Param, Tag, Rest>& params)
williamr@2: {
williamr@2:   typedef typename property_value<bgl_named_params<Param,Tag,Rest>,
williamr@2:                                   vertex_displacement_t>::type D;
williamr@2: 
williamr@2:   detail::fr_force_directed_layout<D>::run
williamr@2:     (g, position, width, height,
williamr@2:      choose_param(get_param(params, attractive_force_t()),
williamr@2:                   square_distance_attractive_force()),
williamr@2:      choose_param(get_param(params, repulsive_force_t()),
williamr@2:                   square_distance_repulsive_force()),
williamr@2:      choose_param(get_param(params, force_pairs_t()),
williamr@2:                   make_grid_force_pairs(width, height, position, g)),
williamr@2:      choose_param(get_param(params, cooling_t()),
williamr@2:                   linear_cooling<Dim>(100)),
williamr@2:      get_param(params, vertex_displacement_t()),
williamr@2:      params);
williamr@2: }
williamr@2: 
williamr@2: template<typename Graph, typename PositionMap, typename Dim>
williamr@2: void
williamr@2: fruchterman_reingold_force_directed_layout(const Graph&    g,
williamr@2:                                            PositionMap     position,
williamr@2:                                            Dim             width,
williamr@2:                                            Dim             height)
williamr@2: {
williamr@2:   fruchterman_reingold_force_directed_layout
williamr@2:     (g, position, width, height,
williamr@2:      attractive_force(square_distance_attractive_force()));
williamr@2: }
williamr@2: 
williamr@2: } // end namespace boost
williamr@2: 
williamr@2: #endif // BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP