Chris@16: // Copyright 2004, 2005 The Trustees of Indiana University.
Chris@16: 
Chris@16: // Distributed under the Boost Software License, Version 1.0.
Chris@16: // (See accompanying file LICENSE_1_0.txt or copy at
Chris@16: // http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: //  Authors: Douglas Gregor
Chris@16: //           Andrew Lumsdaine
Chris@16: #ifndef BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP
Chris@16: #define BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP
Chris@16: 
Chris@16: #include <boost/config/no_tr1/cmath.hpp>
Chris@16: #include <boost/graph/graph_traits.hpp>
Chris@16: #include <boost/graph/named_function_params.hpp>
Chris@16: #include <boost/graph/iteration_macros.hpp>
Chris@16: #include <boost/graph/topology.hpp> // For topology concepts
Chris@16: #include <vector>
Chris@16: #include <list>
Chris@16: #include <algorithm> // for std::min and std::max
Chris@16: #include <numeric> // for std::accumulate
Chris@16: #include <cmath> // for std::sqrt and std::fabs
Chris@16: #include <functional>
Chris@16: 
Chris@16: namespace boost {
Chris@16: 
Chris@16: struct square_distance_attractive_force {
Chris@16:   template<typename Graph, typename T>
Chris@16:   T
Chris@16:   operator()(typename graph_traits<Graph>::edge_descriptor,
Chris@16:              T k,
Chris@16:              T d,
Chris@16:              const Graph&) const
Chris@16:   {
Chris@16:     return d * d / k;
Chris@16:   }
Chris@16: };
Chris@16: 
Chris@16: struct square_distance_repulsive_force {
Chris@16:   template<typename Graph, typename T>
Chris@16:   T
Chris@16:   operator()(typename graph_traits<Graph>::vertex_descriptor,
Chris@16:              typename graph_traits<Graph>::vertex_descriptor,
Chris@16:              T k,
Chris@16:              T d,
Chris@16:              const Graph&) const
Chris@16:   {
Chris@16:     return k * k / d;
Chris@16:   }
Chris@16: };
Chris@16: 
Chris@16: template<typename T>
Chris@16: struct linear_cooling {
Chris@16:   typedef T result_type;
Chris@16: 
Chris@16:   linear_cooling(std::size_t iterations)
Chris@16:     : temp(T(iterations) / T(10)), step(0.1) { }
Chris@16: 
Chris@16:   linear_cooling(std::size_t iterations, T temp)
Chris@16:     : temp(temp), step(temp / T(iterations)) { }
Chris@16: 
Chris@16:   T operator()()
Chris@16:   {
Chris@16:     T old_temp = temp;
Chris@16:     temp -= step;
Chris@16:     if (temp < T(0)) temp = T(0);
Chris@16:     return old_temp;
Chris@16:   }
Chris@16: 
Chris@16:  private:
Chris@16:   T temp;
Chris@16:   T step;
Chris@16: };
Chris@16: 
Chris@16: struct all_force_pairs
Chris@16: {
Chris@16:   template<typename Graph, typename ApplyForce >
Chris@16:   void operator()(const Graph& g, ApplyForce apply_force)
Chris@16:   {
Chris@16:     typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
Chris@16:     vertex_iterator v, end;
Chris@16:     for (boost::tie(v, end) = vertices(g); v != end; ++v) {
Chris@16:       vertex_iterator u = v;
Chris@16:       for (++u; u != end; ++u) {
Chris@16:         apply_force(*u, *v);
Chris@16:         apply_force(*v, *u);
Chris@16:       }
Chris@16:     }
Chris@16:   }
Chris@16: };
Chris@16: 
Chris@16: template<typename Topology, typename PositionMap>
Chris@16: struct grid_force_pairs
Chris@16: {
Chris@16:   typedef typename property_traits<PositionMap>::value_type Point;
Chris@16:   BOOST_STATIC_ASSERT (Point::dimensions == 2);
Chris@16:   typedef typename Topology::point_difference_type point_difference_type;
Chris@16: 
Chris@16:   template<typename Graph>
Chris@16:   explicit
Chris@16:   grid_force_pairs(const Topology& topology,
Chris@16:                    PositionMap position, const Graph& g)
Chris@16:     : topology(topology), position(position)
Chris@16:   {
Chris@16:     two_k = 2. * this->topology.volume(this->topology.extent()) / std::sqrt((double)num_vertices(g));
Chris@16:   }
Chris@16: 
Chris@16:   template<typename Graph, typename ApplyForce >
Chris@16:   void operator()(const Graph& g, ApplyForce apply_force)
Chris@16:   {
Chris@16:     typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
Chris@16:     typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
Chris@16:     typedef std::list<vertex_descriptor> bucket_t;
Chris@16:     typedef std::vector<bucket_t> buckets_t;
Chris@16: 
Chris@16:     std::size_t columns = std::size_t(topology.extent()[0] / two_k + 1.);
Chris@16:     std::size_t rows = std::size_t(topology.extent()[1] / two_k + 1.);
Chris@16:     buckets_t buckets(rows * columns);
Chris@16:     vertex_iterator v, v_end;
Chris@16:     for (boost::tie(v, v_end) = vertices(g); v != v_end; ++v) {
Chris@16:       std::size_t column =
Chris@16:         std::size_t((get(position, *v)[0] + topology.extent()[0] / 2) / two_k);
Chris@16:       std::size_t row    =
Chris@16:         std::size_t((get(position, *v)[1] + topology.extent()[1] / 2) / two_k);
Chris@16: 
Chris@16:       if (column >= columns) column = columns - 1;
Chris@16:       if (row >= rows) row = rows - 1;
Chris@16:       buckets[row * columns + column].push_back(*v);
Chris@16:     }
Chris@16: 
Chris@16:     for (std::size_t row = 0; row < rows; ++row)
Chris@16:       for (std::size_t column = 0; column < columns; ++column) {
Chris@16:         bucket_t& bucket = buckets[row * columns + column];
Chris@16:         typedef typename bucket_t::iterator bucket_iterator;
Chris@16:         for (bucket_iterator u = bucket.begin(); u != bucket.end(); ++u) {
Chris@16:           // Repulse vertices in this bucket
Chris@16:           bucket_iterator v = u;
Chris@16:           for (++v; v != bucket.end(); ++v) {
Chris@16:             apply_force(*u, *v);
Chris@16:             apply_force(*v, *u);
Chris@16:           }
Chris@16: 
Chris@16:           std::size_t adj_start_row = row == 0? 0 : row - 1;
Chris@16:           std::size_t adj_end_row = row == rows - 1? row : row + 1;
Chris@16:           std::size_t adj_start_column = column == 0? 0 : column - 1;
Chris@16:           std::size_t adj_end_column = column == columns - 1? column : column + 1;
Chris@16:           for (std::size_t other_row = adj_start_row; other_row <= adj_end_row;
Chris@16:                ++other_row)
Chris@16:             for (std::size_t other_column = adj_start_column;
Chris@16:                  other_column <= adj_end_column; ++other_column)
Chris@16:               if (other_row != row || other_column != column) {
Chris@16:                 // Repulse vertices in this bucket
Chris@16:                 bucket_t& other_bucket
Chris@16:                   = buckets[other_row * columns + other_column];
Chris@16:                 for (v = other_bucket.begin(); v != other_bucket.end(); ++v) {
Chris@16:                   double dist =
Chris@16:                     topology.distance(get(position, *u), get(position, *v));
Chris@16:                   if (dist < two_k) apply_force(*u, *v);
Chris@16:                 }
Chris@16:               }
Chris@16:         }
Chris@16:       }
Chris@16:   }
Chris@16: 
Chris@16:  private:
Chris@16:   const Topology& topology;
Chris@16:   PositionMap position;
Chris@16:   double two_k;
Chris@16: };
Chris@16: 
Chris@16: template<typename PositionMap, typename Topology, typename Graph>
Chris@16: inline grid_force_pairs<Topology, PositionMap>
Chris@16: make_grid_force_pairs
Chris@16:   (const Topology& topology,
Chris@16:    const PositionMap& position, const Graph& g)
Chris@16: { return grid_force_pairs<Topology, PositionMap>(topology, position, g); }
Chris@16: 
Chris@16: template<typename Graph, typename PositionMap, typename Topology>
Chris@16: void
Chris@16: scale_graph(const Graph& g, PositionMap position, const Topology& topology,
Chris@16:             typename Topology::point_type upper_left, typename Topology::point_type lower_right)
Chris@16: {
Chris@16:   if (num_vertices(g) == 0) return;
Chris@16: 
Chris@16:   typedef typename Topology::point_type Point;
Chris@16:   typedef typename Topology::point_difference_type point_difference_type;
Chris@16: 
Chris@16:   // Find min/max ranges
Chris@16:   Point min_point = get(position, *vertices(g).first), max_point = min_point;
Chris@16:   BGL_FORALL_VERTICES_T(v, g, Graph) {
Chris@16:     min_point = topology.pointwise_min(min_point, get(position, v));
Chris@16:     max_point = topology.pointwise_max(max_point, get(position, v));
Chris@16:   }
Chris@16: 
Chris@16:   Point old_origin = topology.move_position_toward(min_point, 0.5, max_point);
Chris@16:   Point new_origin = topology.move_position_toward(upper_left, 0.5, lower_right);
Chris@16:   point_difference_type old_size = topology.difference(max_point, min_point);
Chris@16:   point_difference_type new_size = topology.difference(lower_right, upper_left);
Chris@16: 
Chris@16:   // Scale to bounding box provided
Chris@16:   BGL_FORALL_VERTICES_T(v, g, Graph) {
Chris@16:     point_difference_type relative_loc = topology.difference(get(position, v), old_origin);
Chris@16:     relative_loc = (relative_loc / old_size) * new_size;
Chris@16:     put(position, v, topology.adjust(new_origin, relative_loc));
Chris@16:   }
Chris@16: }
Chris@16: 
Chris@16: namespace detail {
Chris@16:   template<typename Topology, typename PropMap, typename Vertex>
Chris@16:   void 
Chris@16:   maybe_jitter_point(const Topology& topology,
Chris@16:                      const PropMap& pm, Vertex v,
Chris@16:                      const typename Topology::point_type& p2)
Chris@16:   {
Chris@16:     double too_close = topology.norm(topology.extent()) / 10000.;
Chris@16:     if (topology.distance(get(pm, v), p2) < too_close) {
Chris@16:       put(pm, v, 
Chris@16:           topology.move_position_toward(get(pm, v), 1./200,
Chris@16:                                         topology.random_point()));
Chris@16:     }
Chris@16:   }
Chris@16: 
Chris@16:   template<typename Topology, typename PositionMap, typename DisplacementMap,
Chris@16:            typename RepulsiveForce, typename Graph>
Chris@16:   struct fr_apply_force
Chris@16:   {
Chris@16:     typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
Chris@16:     typedef typename Topology::point_type Point;
Chris@16:     typedef typename Topology::point_difference_type PointDiff;
Chris@16: 
Chris@16:     fr_apply_force(const Topology& topology,
Chris@16:                    const PositionMap& position,
Chris@16:                    const DisplacementMap& displacement,
Chris@16:                    RepulsiveForce repulsive_force, double k, const Graph& g)
Chris@16:       : topology(topology), position(position), displacement(displacement),
Chris@16:         repulsive_force(repulsive_force), k(k), g(g)
Chris@16:     { }
Chris@16: 
Chris@16:     void operator()(vertex_descriptor u, vertex_descriptor v)
Chris@16:     {
Chris@16:       if (u != v) {
Chris@16:         // When the vertices land on top of each other, move the
Chris@16:         // first vertex away from the boundaries.
Chris@16:         maybe_jitter_point(topology, position, u, get(position, v));
Chris@16: 
Chris@16:         double dist = topology.distance(get(position, u), get(position, v));
Chris@16:         typename Topology::point_difference_type dispv = get(displacement, v);
Chris@16:         if (dist == 0.) {
Chris@16:           for (std::size_t i = 0; i < Point::dimensions; ++i) {
Chris@16:             dispv[i] += 0.01;
Chris@16:           }
Chris@16:         } else {
Chris@16:           double fr = repulsive_force(u, v, k, dist, g);
Chris@16:           dispv += (fr / dist) *
Chris@16:                    topology.difference(get(position, v), get(position, u));
Chris@16:         }
Chris@16:         put(displacement, v, dispv);
Chris@16:       }
Chris@16:     }
Chris@16: 
Chris@16:   private:
Chris@16:     const Topology& topology;
Chris@16:     PositionMap position;
Chris@16:     DisplacementMap displacement;
Chris@16:     RepulsiveForce repulsive_force;
Chris@16:     double k;
Chris@16:     const Graph& g;
Chris@16:   };
Chris@16: 
Chris@16: } // end namespace detail
Chris@16: 
Chris@16: template<typename Topology, typename Graph, typename PositionMap, 
Chris@16:          typename AttractiveForce, typename RepulsiveForce,
Chris@16:          typename ForcePairs, typename Cooling, typename DisplacementMap>
Chris@16: void
Chris@16: fruchterman_reingold_force_directed_layout
Chris@16:  (const Graph&    g,
Chris@16:   PositionMap     position,
Chris@16:   const Topology& topology,
Chris@16:   AttractiveForce attractive_force,
Chris@16:   RepulsiveForce  repulsive_force,
Chris@16:   ForcePairs      force_pairs,
Chris@16:   Cooling         cool,
Chris@16:   DisplacementMap displacement)
Chris@16: {
Chris@16:   typedef typename graph_traits<Graph>::vertex_iterator   vertex_iterator;
Chris@16:   typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
Chris@16:   typedef typename graph_traits<Graph>::edge_iterator     edge_iterator;
Chris@16: 
Chris@16:   double volume = topology.volume(topology.extent());
Chris@16: 
Chris@16:   // assume positions are initialized randomly
Chris@16:   double k = pow(volume / num_vertices(g), 1. / (double)(Topology::point_difference_type::dimensions));
Chris@16: 
Chris@16:   detail::fr_apply_force<Topology, PositionMap, DisplacementMap,
Chris@16:                          RepulsiveForce, Graph>
Chris@16:     apply_force(topology, position, displacement, repulsive_force, k, g);
Chris@16: 
Chris@16:   do {
Chris@16:     // Calculate repulsive forces
Chris@16:     vertex_iterator v, v_end;
Chris@16:     for (boost::tie(v, v_end) = vertices(g); v != v_end; ++v)
Chris@16:       put(displacement, *v, typename Topology::point_difference_type());
Chris@16:     force_pairs(g, apply_force);
Chris@16: 
Chris@16:     // Calculate attractive forces
Chris@16:     edge_iterator e, e_end;
Chris@16:     for (boost::tie(e, e_end) = edges(g); e != e_end; ++e) {
Chris@16:       vertex_descriptor v = source(*e, g);
Chris@16:       vertex_descriptor u = target(*e, g);
Chris@16: 
Chris@16:       // When the vertices land on top of each other, move the
Chris@16:       // first vertex away from the boundaries.
Chris@16:       ::boost::detail::maybe_jitter_point(topology, position, u, get(position, v));
Chris@16: 
Chris@16:       typename Topology::point_difference_type delta =
Chris@16:         topology.difference(get(position, v), get(position, u));
Chris@16:       double dist = topology.distance(get(position, u), get(position, v));
Chris@16:       double fa = attractive_force(*e, k, dist, g);
Chris@16: 
Chris@16:       put(displacement, v, get(displacement, v) - (fa / dist) * delta);
Chris@16:       put(displacement, u, get(displacement, u) + (fa / dist) * delta);
Chris@16:     }
Chris@16: 
Chris@16:     if (double temp = cool()) {
Chris@16:       // Update positions
Chris@16:       BGL_FORALL_VERTICES_T (v, g, Graph) {
Chris@16:         BOOST_USING_STD_MIN();
Chris@16:         BOOST_USING_STD_MAX();
Chris@16:         double disp_size = topology.norm(get(displacement, v));
Chris@16:         put(position, v, topology.adjust(get(position, v), get(displacement, v) * (min BOOST_PREVENT_MACRO_SUBSTITUTION (disp_size, temp) / disp_size)));
Chris@16:         put(position, v, topology.bound(get(position, v)));
Chris@16:       }
Chris@16:     } else {
Chris@16:       break;
Chris@16:     }
Chris@16:   } while (true);
Chris@16: }
Chris@16: 
Chris@16: namespace detail {
Chris@16:   template<typename DisplacementMap>
Chris@16:   struct fr_force_directed_layout
Chris@16:   {
Chris@16:     template<typename Topology, typename Graph, typename PositionMap, 
Chris@16:              typename AttractiveForce, typename RepulsiveForce,
Chris@16:              typename ForcePairs, typename Cooling,
Chris@16:              typename Param, typename Tag, typename Rest>
Chris@16:     static void
Chris@16:     run(const Graph&    g,
Chris@16:         PositionMap     position,
Chris@16:         const Topology& topology,
Chris@16:         AttractiveForce attractive_force,
Chris@16:         RepulsiveForce  repulsive_force,
Chris@16:         ForcePairs      force_pairs,
Chris@16:         Cooling         cool,
Chris@16:         DisplacementMap displacement,
Chris@16:         const bgl_named_params<Param, Tag, Rest>&)
Chris@16:     {
Chris@16:       fruchterman_reingold_force_directed_layout
Chris@16:         (g, position, topology, attractive_force, repulsive_force,
Chris@16:          force_pairs, cool, displacement);
Chris@16:     }
Chris@16:   };
Chris@16: 
Chris@16:   template<>
Chris@16:   struct fr_force_directed_layout<param_not_found>
Chris@16:   {
Chris@16:     template<typename Topology, typename Graph, typename PositionMap, 
Chris@16:              typename AttractiveForce, typename RepulsiveForce,
Chris@16:              typename ForcePairs, typename Cooling,
Chris@16:              typename Param, typename Tag, typename Rest>
Chris@16:     static void
Chris@16:     run(const Graph&    g,
Chris@16:         PositionMap     position,
Chris@16:         const Topology& topology,
Chris@16:         AttractiveForce attractive_force,
Chris@16:         RepulsiveForce  repulsive_force,
Chris@16:         ForcePairs      force_pairs,
Chris@16:         Cooling         cool,
Chris@16:         param_not_found,
Chris@16:         const bgl_named_params<Param, Tag, Rest>& params)
Chris@16:     {
Chris@16:       typedef typename Topology::point_difference_type PointDiff;
Chris@16:       std::vector<PointDiff> displacements(num_vertices(g));
Chris@16:       fruchterman_reingold_force_directed_layout
Chris@16:         (g, position, topology, attractive_force, repulsive_force,
Chris@16:          force_pairs, cool,
Chris@16:          make_iterator_property_map
Chris@16:          (displacements.begin(),
Chris@16:           choose_const_pmap(get_param(params, vertex_index), g,
Chris@16:                             vertex_index),
Chris@16:           PointDiff()));
Chris@16:     }
Chris@16:   };
Chris@16: 
Chris@16: } // end namespace detail
Chris@16: 
Chris@16: template<typename Topology, typename Graph, typename PositionMap, typename Param,
Chris@16:          typename Tag, typename Rest>
Chris@16: void
Chris@16: fruchterman_reingold_force_directed_layout
Chris@16:   (const Graph&    g,
Chris@16:    PositionMap     position,
Chris@16:    const Topology& topology,
Chris@16:    const bgl_named_params<Param, Tag, Rest>& params)
Chris@16: {
Chris@16:   typedef typename get_param_type<vertex_displacement_t, bgl_named_params<Param,Tag,Rest> >::type D;
Chris@16: 
Chris@16:   detail::fr_force_directed_layout<D>::run
Chris@16:     (g, position, topology, 
Chris@16:      choose_param(get_param(params, attractive_force_t()),
Chris@16:                   square_distance_attractive_force()),
Chris@16:      choose_param(get_param(params, repulsive_force_t()),
Chris@16:                   square_distance_repulsive_force()),
Chris@16:      choose_param(get_param(params, force_pairs_t()),
Chris@16:                   make_grid_force_pairs(topology, position, g)),
Chris@16:      choose_param(get_param(params, cooling_t()),
Chris@16:                   linear_cooling<double>(100)),
Chris@16:      get_param(params, vertex_displacement_t()),
Chris@16:      params);
Chris@16: }
Chris@16: 
Chris@16: template<typename Topology, typename Graph, typename PositionMap>
Chris@16: void
Chris@16: fruchterman_reingold_force_directed_layout
Chris@16:   (const Graph&    g,
Chris@16:    PositionMap     position,
Chris@16:    const Topology& topology)
Chris@16: {
Chris@16:   fruchterman_reingold_force_directed_layout
Chris@16:     (g, position, topology,
Chris@16:      attractive_force(square_distance_attractive_force()));
Chris@16: }
Chris@16: 
Chris@16: } // end namespace boost
Chris@16: 
Chris@16: #ifdef BOOST_GRAPH_USE_MPI
Chris@16: #  include <boost/graph/distributed/fruchterman_reingold.hpp>
Chris@16: #endif
Chris@16: 
Chris@16: #endif // BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP