boost/graph/fruchterman_reingold.hpp

   1 // Copyright 2004, 2005 The Trustees of Indiana University.
   2
   3 // Distributed under the Boost Software License, Version 1.0.
   4 // (See accompanying file LICENSE_1_0.txt or copy at
   5 // http://www.boost.org/LICENSE_1_0.txt)
   6
   7 //  Authors: Douglas Gregor
   8 //           Andrew Lumsdaine
   9 #ifndef BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP
  10 #define BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP
  11
  12 #include <boost/config/no_tr1/cmath.hpp>
  13 #include <boost/graph/graph_traits.hpp>
  14 #include <boost/graph/named_function_params.hpp>
  15 #include <boost/graph/iteration_macros.hpp>
  16 #include <boost/graph/topology.hpp> // For topology concepts
  17 #include <vector>
  18 #include <list>
  19 #include <algorithm> // for std::min and std::max
  20 #include <numeric> // for std::accumulate
  21 #include <cmath> // for std::sqrt and std::fabs
  22 #include <functional>
  23
  24 namespace boost {
  25
  26   bool vertex_migration = false;
  27
  28 struct square_distance_attractive_force {
  29   template<typename Graph, typename T>
  30   T
  31   operator()(typename graph_traits<Graph>::edge_descriptor,
  32              T k,
  33              T d,
  34              const Graph&) const
  35   {
  36     return d * d / k;
  37   }
  38 };
  39
  40 struct square_distance_repulsive_force {
  41   template<typename Graph, typename T>
  42   T
  43   operator()(typename graph_traits<Graph>::vertex_descriptor,
  44              typename graph_traits<Graph>::vertex_descriptor,
  45              T k,
  46              T d,
  47              const Graph&) const
  48   {
  49     return k * k / d;
  50   }
  51 };
  52
  53 template<typename T>
  54 struct linear_cooling {
  55   typedef T result_type;
  56
  57   linear_cooling(std::size_t iterations)
  58     : temp(T(iterations) / T(10)), step(0.1) { }
  59
  60   linear_cooling(std::size_t iterations, T temp)
  61     : temp(temp), step(temp / T(iterations)) { }
  62
  63   T operator()()
  64   {
  65     T old_temp = temp;
  66     temp -= step;
  67     if (temp < T(0)) temp = T(0);
  68     return old_temp;
  69   }
  70
  71  private:
  72   T temp;
  73   T step;
  74 };
  75
  76 struct all_force_pairs
  77 {
  78   template<typename Graph, typename ApplyForce >
  79   void operator()(const Graph& g, ApplyForce apply_force)
  80   {
  81     typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
  82     vertex_iterator v, end;
  83     for (tie(v, end) = vertices(g); v != end; ++v) {
  84       vertex_iterator u = v;
  85       for (++u; u != end; ++u) {
  86         apply_force(*u, *v);
  87         apply_force(*v, *u);
  88       }
  89     }
  90   }
  91 };
  92
  93 template<typename Topology, typename PositionMap>
  94 struct grid_force_pairs
  95 {
  96   typedef typename property_traits<PositionMap>::value_type Point;
  97   BOOST_STATIC_ASSERT (Point::dimensions == 2);
  98   typedef typename Topology::point_difference_type point_difference_type;
  99
 100   template<typename Graph>
 101   explicit
 102   grid_force_pairs(const Topology& topology,
 103                    const Point& origin, const point_difference_type& extent,
 104                    PositionMap position, const Graph& g)
 105     : topology(topology), extent(extent), origin(origin), position(position)
 106   {
 107     two_k = 2. * this->topology.volume(this->extent) / std::sqrt((double)num_vertices(g));
 108   }
 109
 110   template<typename Graph, typename ApplyForce >
 111   void operator()(const Graph& g, ApplyForce apply_force)
 112   {
 113     typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
 114     typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
 115     typedef std::list<vertex_descriptor> bucket_t;
 116     typedef std::vector<bucket_t> buckets_t;
 117
 118     std::size_t columns = std::size_t(extent[0] / two_k + 1.);
 119     std::size_t rows = std::size_t(extent[1] / two_k + 1.);
 120     buckets_t buckets(rows * columns);
 121     vertex_iterator v, v_end;
 122     for (tie(v, v_end) = vertices(g); v != v_end; ++v) {
 123       std::size_t column = std::size_t((position[*v][0] + extent[0] / 2) / two_k);
 124       std::size_t row    = std::size_t((position[*v][1] + extent[1] / 2) / two_k);
 125
 126       if (column >= columns) column = columns - 1;
 127       if (row >= rows) row = rows - 1;
 128       buckets[row * columns + column].push_back(*v);
 129     }
 130
 131     for (std::size_t row = 0; row < rows; ++row)
 132       for (std::size_t column = 0; column < columns; ++column) {
 133         bucket_t& bucket = buckets[row * columns + column];
 134         typedef typename bucket_t::iterator bucket_iterator;
 135         for (bucket_iterator u = bucket.begin(); u != bucket.end(); ++u) {
 136           // Repulse vertices in this bucket
 137           bucket_iterator v = u;
 138           for (++v; v != bucket.end(); ++v) {
 139             apply_force(*u, *v);
 140             apply_force(*v, *u);
 141           }
 142
 143           std::size_t adj_start_row = row == 0? 0 : row - 1;
 144           std::size_t adj_end_row = row == rows - 1? row : row + 1;
 145           std::size_t adj_start_column = column == 0? 0 : column - 1;
 146           std::size_t adj_end_column = column == columns - 1? column : column + 1;
 147           for (std::size_t other_row = adj_start_row; other_row <= adj_end_row;
 148                ++other_row)
 149             for (std::size_t other_column = adj_start_column;
 150                  other_column <= adj_end_column; ++other_column)
 151               if (other_row != row || other_column != column) {
 152                 // Repulse vertices in this bucket
 153                 bucket_t& other_bucket
 154                   = buckets[other_row * columns + other_column];
 155                 for (v = other_bucket.begin(); v != other_bucket.end(); ++v) {
 156                   double dist = topology.distance(position[*u], position[*v]);
 157                   if (dist < two_k) apply_force(*u, *v);
 158                 }
 159               }
 160         }
 161       }
 162   }
 163
 164  private:
 165   const Topology& topology;
 166   point_difference_type extent;
 167   Point origin;
 168   PositionMap position;
 169   double two_k;
 170 };
 171
 172 template<typename PositionMap, typename Topology, typename Graph>
 173 inline grid_force_pairs<Topology, PositionMap>
 174 make_grid_force_pairs
 175   (const Topology& topology,
 176    typename Topology::point_type const& origin,
 177    typename Topology::point_difference_type const& extent,
 178    const PositionMap& position, const Graph& g)
 179 { return grid_force_pairs<Topology, PositionMap>(topology, origin, extent, position, g); }
 180
 181 template<typename Graph, typename PositionMap, typename Topology>
 182 void
 183 scale_graph(const Graph& g, PositionMap position, const Topology& topology,
 184             typename Topology::point_type upper_left, typename Topology::point_type lower_right)
 185 {
 186   if (num_vertices(g) == 0) return;
 187
 188   typedef typename Topology::point_type Point;
 189   typedef typename Topology::point_difference_type point_difference_type;
 190
 191   // Find min/max ranges
 192   Point min_point = position[*vertices(g).first], max_point = min_point;
 193   BGL_FORALL_VERTICES_T(v, g, Graph) {
 194     min_point = topology.pointwise_min(min_point, position[v]);
 195     max_point = topology.pointwise_max(max_point, position[v]);
 196   }
 197
 198   Point old_origin = topology.move_position_toward(min_point, 0.5, max_point);
 199   Point new_origin = topology.move_position_toward(upper_left, 0.5, lower_right);
 200   point_difference_type old_size = topology.difference(max_point, min_point);
 201   point_difference_type new_size = topology.difference(lower_right, upper_left);
 202
 203   // Scale to bounding box provided
 204   BGL_FORALL_VERTICES_T(v, g, Graph) {
 205     point_difference_type relative_loc = topology.difference(position[v], old_origin);
 206     relative_loc = (relative_loc / old_size) * new_size;
 207     position[v] = topology.adjust(new_origin, relative_loc);
 208   }
 209 }
 210
 211 namespace detail {
 212   template<typename Topology>
 213   void
 214   maybe_jitter_point(const Topology& topology,
 215                      typename Topology::point_type& p1, const typename Topology::point_type& p2,
 216                      typename Topology::point_type origin, typename Topology::point_difference_type extent)
 217   {
 218     double too_close = topology.norm(extent) / 10000.;
 219     if (topology.distance(p1, p2) < too_close) {
 220       p1 = topology.move_position_toward(p1, 1./200, topology.random_point());
 221     }
 222   }
 223
 224   template<typename Topology, typename PositionMap, typename DisplacementMap,
 225            typename RepulsiveForce, typename Graph>
 226   struct fr_apply_force
 227   {
 228     typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
 229     typedef typename Topology::point_type Point;
 230     typedef typename Topology::point_difference_type PointDiff;
 231
 232     fr_apply_force(const Topology& topology,
 233                    const PositionMap& position,
 234                    const DisplacementMap& displacement,
 235                    Point origin, PointDiff extent,
 236                    RepulsiveForce repulsive_force, double k, const Graph& g)
 237       : topology(topology), position(position), displacement(displacement), origin(origin),
 238         extent(extent), repulsive_force(repulsive_force), k(k), g(g)
 239     { }
 240
 241     void operator()(vertex_descriptor u, vertex_descriptor v)
 242     {
 243       if (u != v) {
 244         // When the vertices land on top of each other, move the
 245         // first vertex away from the boundaries.
 246         maybe_jitter_point(topology, position[u], position[v], origin, extent);
 247
 248         double dist = topology.distance(position[u], position[v]);
 249         if (dist == 0.) {
 250           for (std::size_t i = 0; i < Point::dimensions; ++i) {
 251             displacement[v][i] += 0.01;
 252           }
 253         } else {
 254           double fr = repulsive_force(u, v, k, dist, g);
 255           typename Topology::point_difference_type dispv = displacement[v];
 256           dispv += (fr / dist) * topology.difference(position[v], position[u]);
 257         }
 258       }
 259     }
 260
 261   private:
 262     const Topology& topology;
 263     PositionMap position;
 264     DisplacementMap displacement;
 265     Point origin;
 266     PointDiff extent;
 267     RepulsiveForce repulsive_force;
 268     double k;
 269     const Graph& g;
 270   };
 271
 272 } // end namespace detail
 273
 274 template<typename Topology, typename Graph, typename PositionMap,
 275          typename AttractiveForce, typename RepulsiveForce,
 276          typename ForcePairs, typename Cooling, typename DisplacementMap>
 277 void
 278 fruchterman_reingold_force_directed_layout
 279  (const Graph&    g,
 280   PositionMap     position,
 281   const Topology& topology,
 282   typename Topology::point_type const& origin,
 283   typename Topology::point_difference_type const& extent,
 284   AttractiveForce attractive_force,
 285   RepulsiveForce  repulsive_force,
 286   ForcePairs      force_pairs,
 287   Cooling         cool,
 288   DisplacementMap displacement)
 289 {
 290   typedef typename Topology::point_type Point;
 291   typedef typename graph_traits<Graph>::vertex_iterator   vertex_iterator;
 292   typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
 293   typedef typename graph_traits<Graph>::edge_iterator     edge_iterator;
 294
 295   double volume = 1.;
 296   for (std::size_t i = 0; i < Topology::point_difference_type::dimensions; ++i)
 297     volume *= extent[i];
 298
 299   // assume positions are initialized randomly
 300   double k = pow(volume / num_vertices(g), 1. / (double)(Topology::point_difference_type::dimensions));
 301
 302   detail::fr_apply_force<Topology, PositionMap, DisplacementMap,
 303                          RepulsiveForce, Graph>
 304     apply_force(topology, position, displacement, origin, extent, repulsive_force, k, g);
 305
 306   do {
 307     // Calculate repulsive forces
 308     vertex_iterator v, v_end;
 309     for (tie(v, v_end) = vertices(g); v != v_end; ++v)
 310       displacement[*v] = typename Topology::point_difference_type();
 311     force_pairs(g, apply_force);
 312
 313     // Calculate attractive forces
 314     edge_iterator e, e_end;
 315     for (tie(e, e_end) = edges(g); e != e_end; ++e) {
 316       vertex_descriptor v = source(*e, g);
 317       vertex_descriptor u = target(*e, g);
 318
 319       // When the vertices land on top of each other, move the
 320       // first vertex away from the boundaries.
 321       ::boost::detail::maybe_jitter_point(topology, position[u], position[v],
 322                                           origin, extent);
 323
 324       typename Topology::point_difference_type delta = topology.difference(position[v], position[u]);
 325       double dist = topology.distance(position[u], position[v]);
 326       double fa = attractive_force(*e, k, dist, g);
 327
 328       displacement[v] -= (fa / dist) * delta;
 329       displacement[u] += (fa / dist) * delta;
 330     }
 331
 332     if (double temp = cool()) {
 333       // Update positions
 334       BGL_FORALL_VERTICES_T (v, g, Graph) {
 335         BOOST_USING_STD_MIN();
 336         BOOST_USING_STD_MAX();
 337         double disp_size = topology.norm(displacement[v]);
 338         position[v] = topology.adjust(position[v], displacement[v] * (min BOOST_PREVENT_MACRO_SUBSTITUTION (disp_size, temp) / disp_size));
 339         position[v] = topology.bound(position[v]);
 340 #if 0
 341         // CEM HACK: Jitter if we're on the edges
 342         if(position[v].x == 1.0f) // origin.x + extent.x)
 343           position[v].x -= drand48() * .1 * extent.x;
 344         else if(position[v].x == -1.0f) // origin.x)
 345           position[v].x += drand48() * .1 * extent.x;
 346 #endif
 347       }
 348     } else {
 349       break;
 350     }
 351   } while (true);
 352 }
 353
 354 namespace detail {
 355   template<typename DisplacementMap>
 356   struct fr_force_directed_layout
 357   {
 358     template<typename Topology, typename Graph, typename PositionMap,
 359              typename AttractiveForce, typename RepulsiveForce,
 360              typename ForcePairs, typename Cooling,
 361              typename Param, typename Tag, typename Rest>
 362     static void
 363     run(const Graph&    g,
 364         PositionMap     position,
 365         const Topology& topology,
 366         typename Topology::point_type const& origin,
 367         typename Topology::point_difference_type const& extent,
 368         AttractiveForce attractive_force,
 369         RepulsiveForce  repulsive_force,
 370         ForcePairs      force_pairs,
 371         Cooling         cool,
 372         DisplacementMap displacement,
 373         const bgl_named_params<Param, Tag, Rest>&)
 374     {
 375       fruchterman_reingold_force_directed_layout
 376         (g, position, topology, origin, extent, attractive_force, repulsive_force,
 377          force_pairs, cool, displacement);
 378     }
 379   };
 380
 381   template<>
 382   struct fr_force_directed_layout<error_property_not_found>
 383   {
 384     template<typename Topology, typename Graph, typename PositionMap,
 385              typename AttractiveForce, typename RepulsiveForce,
 386              typename ForcePairs, typename Cooling,
 387              typename Param, typename Tag, typename Rest>
 388     static void
 389     run(const Graph&    g,
 390         PositionMap     position,
 391         const Topology& topology,
 392         typename Topology::point_type const& origin,
 393         typename Topology::point_difference_type const& extent,
 394         AttractiveForce attractive_force,
 395         RepulsiveForce  repulsive_force,
 396         ForcePairs      force_pairs,
 397         Cooling         cool,
 398         error_property_not_found,
 399         const bgl_named_params<Param, Tag, Rest>& params)
 400     {
 401       typedef typename Topology::point_difference_type PointDiff;
 402       std::vector<PointDiff> displacements(num_vertices(g));
 403       fruchterman_reingold_force_directed_layout
 404         (g, position, topology, origin, extent, attractive_force, repulsive_force,
 405          force_pairs, cool,
 406          make_iterator_property_map
 407          (displacements.begin(),
 408           choose_const_pmap(get_param(params, vertex_index), g,
 409                             vertex_index),
 410           PointDiff()));
 411     }
 412   };
 413
 414 } // end namespace detail
 415
 416 template<typename Topology, typename Graph, typename PositionMap, typename Param,
 417          typename Tag, typename Rest>
 418 void
 419 fruchterman_reingold_force_directed_layout
 420   (const Graph&    g,
 421    PositionMap     position,
 422    const Topology& topology,
 423    typename Topology::point_type const& origin,
 424    typename Topology::point_difference_type const& extent,
 425    const bgl_named_params<Param, Tag, Rest>& params)
 426 {
 427   typedef typename property_value<bgl_named_params<Param,Tag,Rest>,
 428                                   vertex_displacement_t>::type D;
 429
 430   detail::fr_force_directed_layout<D>::run
 431     (g, position, topology, origin, extent,
 432      choose_param(get_param(params, attractive_force_t()),
 433                   square_distance_attractive_force()),
 434      choose_param(get_param(params, repulsive_force_t()),
 435                   square_distance_repulsive_force()),
 436      choose_param(get_param(params, force_pairs_t()),
 437                   make_grid_force_pairs(topology, origin, extent, position, g)),
 438      choose_param(get_param(params, cooling_t()),
 439                   linear_cooling<double>(100)),
 440      get_param(params, vertex_displacement_t()),
 441      params);
 442 }
 443
 444 template<typename Topology, typename Graph, typename PositionMap>
 445 void
 446 fruchterman_reingold_force_directed_layout
 447   (const Graph&    g,
 448    PositionMap     position,
 449    const Topology& topology,
 450    typename Topology::point_type const& origin,
 451    typename Topology::point_difference_type const& extent)
 452 {
 453   fruchterman_reingold_force_directed_layout
 454     (g, position, topology, origin, extent,
 455      attractive_force(square_distance_attractive_force()));
 456 }
 457
 458 } // end namespace boost
 459
 460 #ifdef BOOST_GRAPH_USE_MPI
 461 #  include <boost/graph/distributed/fruchterman_reingold.hpp>
 462 #endif
 463
 464 #endif // BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP