boost/graph/fruchterman_reingold.hpp
// Copyright 2004, 2005 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 <boost/config/no_tr1/cmath.hpp>
#include <boost/graph/graph_traits.hpp>
#include <boost/graph/named_function_params.hpp>
#include <boost/graph/iteration_macros.hpp>
#include <boost/graph/topology.hpp> // For topology concepts
#include <boost/graph/detail/mpi_include.hpp>
#include <vector>
#include <list>
#include <algorithm> // for std::min and std::max
#include <numeric> // for std::accumulate
#include <cmath> // for std::sqrt and std::fabs
#include <functional>
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 (boost::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 Topology, typename PositionMap>
struct grid_force_pairs
{
typedef typename property_traits<PositionMap>::value_type Point;
BOOST_STATIC_ASSERT (Point::dimensions == 2);
typedef typename Topology::point_difference_type point_difference_type;
template<typename Graph>
explicit
grid_force_pairs(const Topology& topology,
PositionMap position, const Graph& g)
: topology(topology), position(position)
{
two_k = 2. * this->topology.volume(this->topology.extent()) / std::sqrt((double)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(topology.extent()[0] / two_k + 1.);
std::size_t rows = std::size_t(topology.extent()[1] / two_k + 1.);
buckets_t buckets(rows * columns);
vertex_iterator v, v_end;
for (boost::tie(v, v_end) = vertices(g); v != v_end; ++v) {
std::size_t column =
std::size_t((get(position, *v)[0] + topology.extent()[0] / 2) / two_k);
std::size_t row =
std::size_t((get(position, *v)[1] + topology.extent()[1] / 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) {
double dist =
topology.distance(get(position, *u), get(position, *v));
if (dist < two_k) apply_force(*u, *v);
}
}
}
}
}
private:
const Topology& topology;
PositionMap position;
double two_k;
};
template<typename PositionMap, typename Topology, typename Graph>
inline grid_force_pairs<Topology, PositionMap>
make_grid_force_pairs
(const Topology& topology,
const PositionMap& position, const Graph& g)
{ return grid_force_pairs<Topology, PositionMap>(topology, position, g); }
template<typename Graph, typename PositionMap, typename Topology>
void
scale_graph(const Graph& g, PositionMap position, const Topology& topology,
typename Topology::point_type upper_left, typename Topology::point_type lower_right)
{
if (num_vertices(g) == 0) return;
typedef typename Topology::point_type Point;
typedef typename Topology::point_difference_type point_difference_type;
// Find min/max ranges
Point min_point = get(position, *vertices(g).first), max_point = min_point;
BGL_FORALL_VERTICES_T(v, g, Graph) {
min_point = topology.pointwise_min(min_point, get(position, v));
max_point = topology.pointwise_max(max_point, get(position, v));
}
Point old_origin = topology.move_position_toward(min_point, 0.5, max_point);
Point new_origin = topology.move_position_toward(upper_left, 0.5, lower_right);
point_difference_type old_size = topology.difference(max_point, min_point);
point_difference_type new_size = topology.difference(lower_right, upper_left);
// Scale to bounding box provided
BGL_FORALL_VERTICES_T(v, g, Graph) {
point_difference_type relative_loc = topology.difference(get(position, v), old_origin);
relative_loc = (relative_loc / old_size) * new_size;
put(position, v, topology.adjust(new_origin, relative_loc));
}
}
namespace detail {
template<typename Topology, typename PropMap, typename Vertex>
void
maybe_jitter_point(const Topology& topology,
const PropMap& pm, Vertex v,
const typename Topology::point_type& p2)
{
double too_close = topology.norm(topology.extent()) / 10000.;
if (topology.distance(get(pm, v), p2) < too_close) {
put(pm, v,
topology.move_position_toward(get(pm, v), 1./200,
topology.random_point()));
}
}
template<typename Topology, typename PositionMap, typename DisplacementMap,
typename RepulsiveForce, typename Graph>
struct fr_apply_force
{
typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
typedef typename Topology::point_type Point;
typedef typename Topology::point_difference_type PointDiff;
fr_apply_force(const Topology& topology,
const PositionMap& position,
const DisplacementMap& displacement,
RepulsiveForce repulsive_force, double k, const Graph& g)
: topology(topology), position(position), displacement(displacement),
repulsive_force(repulsive_force), k(k), g(g)
{ }
void operator()(vertex_descriptor u, vertex_descriptor v)
{
if (u != v) {
// When the vertices land on top of each other, move the
// first vertex away from the boundaries.
maybe_jitter_point(topology, position, u, get(position, v));
double dist = topology.distance(get(position, u), get(position, v));
typename Topology::point_difference_type dispv = get(displacement, v);
if (dist == 0.) {
for (std::size_t i = 0; i < Point::dimensions; ++i) {
dispv[i] += 0.01;
}
} else {
double fr = repulsive_force(u, v, k, dist, g);
dispv += (fr / dist) *
topology.difference(get(position, v), get(position, u));
}
put(displacement, v, dispv);
}
}
private:
const Topology& topology;
PositionMap position;
DisplacementMap displacement;
RepulsiveForce repulsive_force;
double k;
const Graph& g;
};
} // end namespace detail
template<typename Topology, typename Graph, typename PositionMap,
typename AttractiveForce, typename RepulsiveForce,
typename ForcePairs, typename Cooling, typename DisplacementMap>
void
fruchterman_reingold_force_directed_layout
(const Graph& g,
PositionMap position,
const Topology& topology,
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;
double volume = topology.volume(topology.extent());
// assume positions are initialized randomly
double k = pow(volume / num_vertices(g), 1. / (double)(Topology::point_difference_type::dimensions));
detail::fr_apply_force<Topology, PositionMap, DisplacementMap,
RepulsiveForce, Graph>
apply_force(topology, position, displacement, repulsive_force, k, g);
do {
// Calculate repulsive forces
vertex_iterator v, v_end;
for (boost::tie(v, v_end) = vertices(g); v != v_end; ++v)
put(displacement, *v, typename Topology::point_difference_type());
force_pairs(g, apply_force);
// Calculate attractive forces
edge_iterator e, e_end;
for (boost::tie(e, e_end) = edges(g); e != e_end; ++e) {
vertex_descriptor v = source(*e, g);
vertex_descriptor u = target(*e, g);
// When the vertices land on top of each other, move the
// first vertex away from the boundaries.
::boost::detail::maybe_jitter_point(topology, position, u, get(position, v));
typename Topology::point_difference_type delta =
topology.difference(get(position, v), get(position, u));
double dist = topology.distance(get(position, u), get(position, v));
double fa = attractive_force(*e, k, dist, g);
put(displacement, v, get(displacement, v) - (fa / dist) * delta);
put(displacement, u, get(displacement, u) + (fa / dist) * delta);
}
if (double temp = cool()) {
// Update positions
BGL_FORALL_VERTICES_T (v, g, Graph) {
BOOST_USING_STD_MIN();
BOOST_USING_STD_MAX();
double disp_size = topology.norm(get(displacement, v));
put(position, v, topology.adjust(get(position, v), get(displacement, v) * (min BOOST_PREVENT_MACRO_SUBSTITUTION (disp_size, temp) / disp_size)));
put(position, v, topology.bound(get(position, v)));
}
} else {
break;
}
} while (true);
}
namespace detail {
template<typename DisplacementMap>
struct fr_force_directed_layout
{
template<typename Topology, typename Graph, typename PositionMap,
typename AttractiveForce, typename RepulsiveForce,
typename ForcePairs, typename Cooling,
typename Param, typename Tag, typename Rest>
static void
run(const Graph& g,
PositionMap position,
const Topology& topology,
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, topology, attractive_force, repulsive_force,
force_pairs, cool, displacement);
}
};
template<>
struct fr_force_directed_layout<param_not_found>
{
template<typename Topology, typename Graph, typename PositionMap,
typename AttractiveForce, typename RepulsiveForce,
typename ForcePairs, typename Cooling,
typename Param, typename Tag, typename Rest>
static void
run(const Graph& g,
PositionMap position,
const Topology& topology,
AttractiveForce attractive_force,
RepulsiveForce repulsive_force,
ForcePairs force_pairs,
Cooling cool,
param_not_found,
const bgl_named_params<Param, Tag, Rest>& params)
{
typedef typename Topology::point_difference_type PointDiff;
std::vector<PointDiff> displacements(num_vertices(g));
fruchterman_reingold_force_directed_layout
(g, position, topology, attractive_force, repulsive_force,
force_pairs, cool,
make_iterator_property_map
(displacements.begin(),
choose_const_pmap(get_param(params, vertex_index), g,
vertex_index),
PointDiff()));
}
};
} // end namespace detail
template<typename Topology, typename Graph, typename PositionMap, typename Param,
typename Tag, typename Rest>
void
fruchterman_reingold_force_directed_layout
(const Graph& g,
PositionMap position,
const Topology& topology,
const bgl_named_params<Param, Tag, Rest>& params)
{
typedef typename get_param_type<vertex_displacement_t, bgl_named_params<Param,Tag,Rest> >::type D;
detail::fr_force_directed_layout<D>::run
(g, position, topology,
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(topology, position, g)),
choose_param(get_param(params, cooling_t()),
linear_cooling<double>(100)),
get_param(params, vertex_displacement_t()),
params);
}
template<typename Topology, typename Graph, typename PositionMap>
void
fruchterman_reingold_force_directed_layout
(const Graph& g,
PositionMap position,
const Topology& topology)
{
fruchterman_reingold_force_directed_layout
(g, position, topology,
attractive_force(square_distance_attractive_force()));
}
} // end namespace boost
#include BOOST_GRAPH_MPI_INCLUDE(<boost/graph/distributed/fruchterman_reingold.hpp>)
#endif // BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP