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