| // Copyright 2004 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_KAMADA_KAWAI_SPRING_LAYOUT_HPP |
| #define BOOST_GRAPH_KAMADA_KAWAI_SPRING_LAYOUT_HPP |
| |
| #include <boost/graph/graph_traits.hpp> |
| #include <boost/graph/topology.hpp> |
| #include <boost/graph/iteration_macros.hpp> |
| #include <boost/graph/johnson_all_pairs_shortest.hpp> |
| #include <boost/type_traits/is_convertible.hpp> |
| #include <utility> |
| #include <iterator> |
| #include <vector> |
| #include <boost/limits.hpp> |
| #include <boost/config/no_tr1/cmath.hpp> |
| |
| namespace boost { |
| namespace detail { namespace graph { |
| /** |
| * Denotes an edge or display area side length used to scale a |
| * Kamada-Kawai drawing. |
| */ |
| template<bool Edge, typename T> |
| struct edge_or_side |
| { |
| explicit edge_or_side(T value) : value(value) {} |
| |
| T value; |
| }; |
| |
| /** |
| * Compute the edge length from an edge length. This is trivial. |
| */ |
| template<typename Graph, typename DistanceMap, typename IndexMap, |
| typename T> |
| T compute_edge_length(const Graph&, DistanceMap, IndexMap, |
| edge_or_side<true, T> length) |
| { return length.value; } |
| |
| /** |
| * Compute the edge length based on the display area side |
| length. We do this by dividing the side length by the largest |
| shortest distance between any two vertices in the graph. |
| */ |
| template<typename Graph, typename DistanceMap, typename IndexMap, |
| typename T> |
| T |
| compute_edge_length(const Graph& g, DistanceMap distance, IndexMap index, |
| edge_or_side<false, T> length) |
| { |
| T result(0); |
| |
| typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator; |
| |
| for (vertex_iterator ui = vertices(g).first, end = vertices(g).second; |
| ui != end; ++ui) { |
| vertex_iterator vi = ui; |
| for (++vi; vi != end; ++vi) { |
| T dij = distance[get(index, *ui)][get(index, *vi)]; |
| if (dij > result) result = dij; |
| } |
| } |
| return length.value / result; |
| } |
| |
| /** |
| * Dense linear solver for fixed-size matrices. |
| */ |
| template <std::size_t Size> |
| struct linear_solver { |
| // Indices in mat are (row, column) |
| // template <typename Vec> |
| // static Vec solve(double mat[Size][Size], Vec rhs); |
| }; |
| |
| template <> |
| struct linear_solver<1> { |
| template <typename Vec> |
| static Vec solve(double mat[1][1], Vec rhs) { |
| return rhs / mat[0][0]; |
| } |
| }; |
| |
| // These are from http://en.wikipedia.org/wiki/Cramer%27s_rule |
| template <> |
| struct linear_solver<2> { |
| template <typename Vec> |
| static Vec solve(double mat[2][2], Vec rhs) { |
| double denom = mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1]; |
| double x_num = rhs[0] * mat[1][1] - rhs[1] * mat[0][1]; |
| double y_num = mat[0][0] * rhs[1] - mat[1][0] * rhs[0] ; |
| Vec result; |
| result[0] = x_num / denom; |
| result[1] = y_num / denom; |
| return result; |
| } |
| }; |
| |
| template <> |
| struct linear_solver<3> { |
| template <typename Vec> |
| static Vec solve(double mat[2][2], Vec rhs) { |
| double denom = mat[0][0] * (mat[1][1] * mat[2][2] - mat[2][1] * mat[1][2]) |
| - mat[1][0] * (mat[0][1] * mat[2][2] - mat[2][1] * mat[0][2]) |
| + mat[2][0] * (mat[0][1] * mat[1][2] - mat[1][1] * mat[0][2]); |
| double x_num = rhs[0] * (mat[1][1] * mat[2][2] - mat[2][1] * mat[1][2]) |
| - rhs[1] * (mat[0][1] * mat[2][2] - mat[2][1] * mat[0][2]) |
| + rhs[2] * (mat[0][1] * mat[1][2] - mat[1][1] * mat[0][2]); |
| double y_num = mat[0][0] * (rhs[1] * mat[2][2] - rhs[2] * mat[1][2]) |
| - mat[1][0] * (rhs[0] * mat[2][2] - rhs[2] * mat[0][2]) |
| + mat[2][0] * (rhs[0] * mat[1][2] - rhs[1] * mat[0][2]); |
| double z_num = mat[0][0] * (mat[1][1] * rhs[2] - mat[2][1] * rhs[1] ) |
| - mat[1][0] * (mat[0][1] * rhs[2] - mat[2][1] * rhs[0] ) |
| + mat[2][0] * (mat[0][1] * rhs[1] - mat[1][1] * rhs[0] ); |
| Vec result; |
| result[0] = x_num / denom; |
| result[1] = y_num / denom; |
| result[2] = z_num / denom; |
| return result; |
| } |
| }; |
| |
| /** |
| * Implementation of the Kamada-Kawai spring layout algorithm. |
| */ |
| template<typename Topology, typename Graph, typename PositionMap, typename WeightMap, |
| typename EdgeOrSideLength, typename Done, |
| typename VertexIndexMap, typename DistanceMatrix, |
| typename SpringStrengthMatrix, typename PartialDerivativeMap> |
| struct kamada_kawai_spring_layout_impl |
| { |
| typedef typename property_traits<WeightMap>::value_type weight_type; |
| typedef typename Topology::point_type Point; |
| typedef typename Topology::point_difference_type point_difference_type; |
| typedef point_difference_type deriv_type; |
| typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator; |
| typedef typename graph_traits<Graph>::vertex_descriptor |
| vertex_descriptor; |
| |
| kamada_kawai_spring_layout_impl( |
| const Topology& topology, |
| const Graph& g, |
| PositionMap position, |
| WeightMap weight, |
| EdgeOrSideLength edge_or_side_length, |
| Done done, |
| weight_type spring_constant, |
| VertexIndexMap index, |
| DistanceMatrix distance, |
| SpringStrengthMatrix spring_strength, |
| PartialDerivativeMap partial_derivatives) |
| : topology(topology), g(g), position(position), weight(weight), |
| edge_or_side_length(edge_or_side_length), done(done), |
| spring_constant(spring_constant), index(index), distance(distance), |
| spring_strength(spring_strength), |
| partial_derivatives(partial_derivatives) {} |
| |
| // Compute contribution of vertex i to the first partial |
| // derivatives (dE/dx_m, dE/dy_m) (for vertex m) |
| deriv_type |
| compute_partial_derivative(vertex_descriptor m, vertex_descriptor i) |
| { |
| #ifndef BOOST_NO_STDC_NAMESPACE |
| using std::sqrt; |
| #endif // BOOST_NO_STDC_NAMESPACE |
| |
| deriv_type result; |
| if (i != m) { |
| point_difference_type diff = topology.difference(position[m], position[i]); |
| weight_type dist = topology.norm(diff); |
| result = spring_strength[get(index, m)][get(index, i)] |
| * (diff - distance[get(index, m)][get(index, i)]/dist*diff); |
| } |
| |
| return result; |
| } |
| |
| // Compute partial derivatives dE/dx_m and dE/dy_m |
| deriv_type |
| compute_partial_derivatives(vertex_descriptor m) |
| { |
| #ifndef BOOST_NO_STDC_NAMESPACE |
| using std::sqrt; |
| #endif // BOOST_NO_STDC_NAMESPACE |
| |
| deriv_type result; |
| |
| // TBD: looks like an accumulate to me |
| BGL_FORALL_VERTICES_T(i, g, Graph) { |
| deriv_type deriv = compute_partial_derivative(m, i); |
| result += deriv; |
| } |
| |
| return result; |
| } |
| |
| // The actual Kamada-Kawai spring layout algorithm implementation |
| bool run() |
| { |
| #ifndef BOOST_NO_STDC_NAMESPACE |
| using std::sqrt; |
| #endif // BOOST_NO_STDC_NAMESPACE |
| |
| // Compute d_{ij} and place it in the distance matrix |
| if (!johnson_all_pairs_shortest_paths(g, distance, index, weight, |
| weight_type(0))) |
| return false; |
| |
| // Compute L based on side length (if needed), or retrieve L |
| weight_type edge_length = |
| detail::graph::compute_edge_length(g, distance, index, |
| edge_or_side_length); |
| |
| // Compute l_{ij} and k_{ij} |
| const weight_type K = spring_constant; |
| vertex_iterator ui, end = vertices(g).second; |
| for (ui = vertices(g).first; ui != end; ++ui) { |
| vertex_iterator vi = ui; |
| for (++vi; vi != end; ++vi) { |
| weight_type dij = distance[get(index, *ui)][get(index, *vi)]; |
| if (dij == (std::numeric_limits<weight_type>::max)()) |
| return false; |
| distance[get(index, *ui)][get(index, *vi)] = edge_length * dij; |
| distance[get(index, *vi)][get(index, *ui)] = edge_length * dij; |
| spring_strength[get(index, *ui)][get(index, *vi)] = K/(dij*dij); |
| spring_strength[get(index, *vi)][get(index, *ui)] = K/(dij*dij); |
| } |
| } |
| |
| // Compute Delta_i and find max |
| vertex_descriptor p = *vertices(g).first; |
| weight_type delta_p(0); |
| |
| for (ui = vertices(g).first; ui != end; ++ui) { |
| deriv_type deriv = compute_partial_derivatives(*ui); |
| put(partial_derivatives, *ui, deriv); |
| |
| weight_type delta = topology.norm(deriv); |
| |
| if (delta > delta_p) { |
| p = *ui; |
| delta_p = delta; |
| } |
| } |
| |
| while (!done(delta_p, p, g, true)) { |
| // The contribution p makes to the partial derivatives of |
| // each vertex. Computing this (at O(n) cost) allows us to |
| // update the delta_i values in O(n) time instead of O(n^2) |
| // time. |
| std::vector<deriv_type> p_partials(num_vertices(g)); |
| for (ui = vertices(g).first; ui != end; ++ui) { |
| vertex_descriptor i = *ui; |
| p_partials[get(index, i)] = compute_partial_derivative(i, p); |
| } |
| |
| do { |
| // Compute the elements of the Jacobian |
| // From |
| // http://www.cs.panam.edu/~rfowler/papers/1994_kumar_fowler_A_Spring_UTPACSTR.pdf |
| // with the bugs fixed in the off-diagonal case |
| weight_type dE_d_d[Point::dimensions][Point::dimensions]; |
| for (std::size_t i = 0; i < Point::dimensions; ++i) |
| for (std::size_t j = 0; j < Point::dimensions; ++j) |
| dE_d_d[i][j] = 0.; |
| for (ui = vertices(g).first; ui != end; ++ui) { |
| vertex_descriptor i = *ui; |
| if (i != p) { |
| point_difference_type diff = topology.difference(position[p], position[i]); |
| weight_type dist = topology.norm(diff); |
| weight_type dist_squared = dist * dist; |
| weight_type inv_dist_cubed = 1. / (dist_squared * dist); |
| weight_type k_mi = spring_strength[get(index,p)][get(index,i)]; |
| weight_type l_mi = distance[get(index, p)][get(index, i)]; |
| for (std::size_t i = 0; i < Point::dimensions; ++i) { |
| for (std::size_t j = 0; j < Point::dimensions; ++j) { |
| if (i == j) { |
| dE_d_d[i][i] += k_mi * (1 + (l_mi * (diff[i] * diff[i] - dist_squared) * inv_dist_cubed)); |
| } else { |
| dE_d_d[i][j] += k_mi * l_mi * diff[i] * diff[j] * inv_dist_cubed; |
| } |
| } |
| } |
| } |
| } |
| |
| deriv_type dE_d = get(partial_derivatives, p); |
| |
| // Solve dE_d_d * delta = dE_d to get delta |
| point_difference_type delta = -linear_solver<Point::dimensions>::solve(dE_d_d, dE_d); |
| |
| // Move p by delta |
| position[p] = topology.adjust(position[p], delta); |
| |
| // Recompute partial derivatives and delta_p |
| deriv_type deriv = compute_partial_derivatives(p); |
| put(partial_derivatives, p, deriv); |
| |
| delta_p = topology.norm(deriv); |
| } while (!done(delta_p, p, g, false)); |
| |
| // Select new p by updating each partial derivative and delta |
| vertex_descriptor old_p = p; |
| for (ui = vertices(g).first; ui != end; ++ui) { |
| deriv_type old_deriv_p = p_partials[get(index, *ui)]; |
| deriv_type old_p_partial = |
| compute_partial_derivative(*ui, old_p); |
| deriv_type deriv = get(partial_derivatives, *ui); |
| |
| deriv += old_p_partial - old_deriv_p; |
| |
| put(partial_derivatives, *ui, deriv); |
| weight_type delta = topology.norm(deriv); |
| |
| if (delta > delta_p) { |
| p = *ui; |
| delta_p = delta; |
| } |
| } |
| } |
| |
| return true; |
| } |
| |
| const Topology& topology; |
| const Graph& g; |
| PositionMap position; |
| WeightMap weight; |
| EdgeOrSideLength edge_or_side_length; |
| Done done; |
| weight_type spring_constant; |
| VertexIndexMap index; |
| DistanceMatrix distance; |
| SpringStrengthMatrix spring_strength; |
| PartialDerivativeMap partial_derivatives; |
| }; |
| } } // end namespace detail::graph |
| |
| /// States that the given quantity is an edge length. |
| template<typename T> |
| inline detail::graph::edge_or_side<true, T> |
| edge_length(T x) |
| { return detail::graph::edge_or_side<true, T>(x); } |
| |
| /// States that the given quantity is a display area side length. |
| template<typename T> |
| inline detail::graph::edge_or_side<false, T> |
| side_length(T x) |
| { return detail::graph::edge_or_side<false, T>(x); } |
| |
| /** |
| * \brief Determines when to terminate layout of a particular graph based |
| * on a given relative tolerance. |
| */ |
| template<typename T = double> |
| struct layout_tolerance |
| { |
| layout_tolerance(const T& tolerance = T(0.001)) |
| : tolerance(tolerance), last_energy((std::numeric_limits<T>::max)()), |
| last_local_energy((std::numeric_limits<T>::max)()) { } |
| |
| template<typename Graph> |
| bool |
| operator()(T delta_p, |
| typename boost::graph_traits<Graph>::vertex_descriptor p, |
| const Graph& g, |
| bool global) |
| { |
| if (global) { |
| if (last_energy == (std::numeric_limits<T>::max)()) { |
| last_energy = delta_p; |
| return false; |
| } |
| |
| T diff = last_energy - delta_p; |
| if (diff < T(0)) diff = -diff; |
| bool done = (delta_p == T(0) || diff / last_energy < tolerance); |
| last_energy = delta_p; |
| return done; |
| } else { |
| if (last_local_energy == (std::numeric_limits<T>::max)()) { |
| last_local_energy = delta_p; |
| return delta_p == T(0); |
| } |
| |
| T diff = last_local_energy - delta_p; |
| bool done = (delta_p == T(0) || (diff / last_local_energy) < tolerance); |
| last_local_energy = delta_p; |
| return done; |
| } |
| } |
| |
| private: |
| T tolerance; |
| T last_energy; |
| T last_local_energy; |
| }; |
| |
| /** \brief Kamada-Kawai spring layout for undirected graphs. |
| * |
| * This algorithm performs graph layout (in two dimensions) for |
| * connected, undirected graphs. It operates by relating the layout |
| * of graphs to a dynamic spring system and minimizing the energy |
| * within that system. The strength of a spring between two vertices |
| * is inversely proportional to the square of the shortest distance |
| * (in graph terms) between those two vertices. Essentially, |
| * vertices that are closer in the graph-theoretic sense (i.e., by |
| * following edges) will have stronger springs and will therefore be |
| * placed closer together. |
| * |
| * Prior to invoking this algorithm, it is recommended that the |
| * vertices be placed along the vertices of a regular n-sided |
| * polygon. |
| * |
| * \param g (IN) must be a model of Vertex List Graph, Edge List |
| * Graph, and Incidence Graph and must be undirected. |
| * |
| * \param position (OUT) must be a model of Lvalue Property Map, |
| * where the value type is a class containing fields @c x and @c y |
| * that will be set to the @c x and @c y coordinates of each vertex. |
| * |
| * \param weight (IN) must be a model of Readable Property Map, |
| * which provides the weight of each edge in the graph @p g. |
| * |
| * \param topology (IN) must be a topology object (see topology.hpp), |
| * which provides operations on points and differences between them. |
| * |
| * \param edge_or_side_length (IN) provides either the unit length |
| * @c e of an edge in the layout or the length of a side @c s of the |
| * display area, and must be either @c boost::edge_length(e) or @c |
| * boost::side_length(s), respectively. |
| * |
| * \param done (IN) is a 4-argument function object that is passed |
| * the current value of delta_p (i.e., the energy of vertex @p p), |
| * the vertex @p p, the graph @p g, and a boolean flag indicating |
| * whether @p delta_p is the maximum energy in the system (when @c |
| * true) or the energy of the vertex being moved. Defaults to @c |
| * layout_tolerance instantiated over the value type of the weight |
| * map. |
| * |
| * \param spring_constant (IN) is the constant multiplied by each |
| * spring's strength. Larger values create systems with more energy |
| * that can take longer to stabilize; smaller values create systems |
| * with less energy that stabilize quickly but do not necessarily |
| * result in pleasing layouts. The default value is 1. |
| * |
| * \param index (IN) is a mapping from vertices to index values |
| * between 0 and @c num_vertices(g). The default is @c |
| * get(vertex_index,g). |
| * |
| * \param distance (UTIL/OUT) will be used to store the distance |
| * from every vertex to every other vertex, which is computed in the |
| * first stages of the algorithm. This value's type must be a model |
| * of BasicMatrix with value type equal to the value type of the |
| * weight map. The default is a a vector of vectors. |
| * |
| * \param spring_strength (UTIL/OUT) will be used to store the |
| * strength of the spring between every pair of vertices. This |
| * value's type must be a model of BasicMatrix with value type equal |
| * to the value type of the weight map. The default is a a vector of |
| * vectors. |
| * |
| * \param partial_derivatives (UTIL) will be used to store the |
| * partial derivates of each vertex with respect to the @c x and @c |
| * y coordinates. This must be a Read/Write Property Map whose value |
| * type is a pair with both types equivalent to the value type of |
| * the weight map. The default is an iterator property map. |
| * |
| * \returns @c true if layout was successful or @c false if a |
| * negative weight cycle was detected. |
| */ |
| template<typename Topology, typename Graph, typename PositionMap, typename WeightMap, |
| typename T, bool EdgeOrSideLength, typename Done, |
| typename VertexIndexMap, typename DistanceMatrix, |
| typename SpringStrengthMatrix, typename PartialDerivativeMap> |
| bool |
| kamada_kawai_spring_layout( |
| const Graph& g, |
| PositionMap position, |
| WeightMap weight, |
| const Topology& topology, |
| detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length, |
| Done done, |
| typename property_traits<WeightMap>::value_type spring_constant, |
| VertexIndexMap index, |
| DistanceMatrix distance, |
| SpringStrengthMatrix spring_strength, |
| PartialDerivativeMap partial_derivatives) |
| { |
| BOOST_STATIC_ASSERT((is_convertible< |
| typename graph_traits<Graph>::directed_category*, |
| undirected_tag* |
| >::value)); |
| |
| detail::graph::kamada_kawai_spring_layout_impl< |
| Topology, Graph, PositionMap, WeightMap, |
| detail::graph::edge_or_side<EdgeOrSideLength, T>, Done, VertexIndexMap, |
| DistanceMatrix, SpringStrengthMatrix, PartialDerivativeMap> |
| alg(topology, g, position, weight, edge_or_side_length, done, spring_constant, |
| index, distance, spring_strength, partial_derivatives); |
| return alg.run(); |
| } |
| |
| /** |
| * \overload |
| */ |
| template<typename Topology, typename Graph, typename PositionMap, typename WeightMap, |
| typename T, bool EdgeOrSideLength, typename Done, |
| typename VertexIndexMap> |
| bool |
| kamada_kawai_spring_layout( |
| const Graph& g, |
| PositionMap position, |
| WeightMap weight, |
| const Topology& topology, |
| detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length, |
| Done done, |
| typename property_traits<WeightMap>::value_type spring_constant, |
| VertexIndexMap index) |
| { |
| typedef typename property_traits<WeightMap>::value_type weight_type; |
| |
| typename graph_traits<Graph>::vertices_size_type n = num_vertices(g); |
| typedef std::vector<weight_type> weight_vec; |
| |
| std::vector<weight_vec> distance(n, weight_vec(n)); |
| std::vector<weight_vec> spring_strength(n, weight_vec(n)); |
| std::vector<typename Topology::point_difference_type> partial_derivatives(n); |
| |
| return |
| kamada_kawai_spring_layout( |
| g, position, weight, topology, edge_or_side_length, done, spring_constant, index, |
| distance.begin(), |
| spring_strength.begin(), |
| make_iterator_property_map(partial_derivatives.begin(), index, |
| typename Topology::point_difference_type())); |
| } |
| |
| /** |
| * \overload |
| */ |
| template<typename Topology, typename Graph, typename PositionMap, typename WeightMap, |
| typename T, bool EdgeOrSideLength, typename Done> |
| bool |
| kamada_kawai_spring_layout( |
| const Graph& g, |
| PositionMap position, |
| WeightMap weight, |
| const Topology& topology, |
| detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length, |
| Done done, |
| typename property_traits<WeightMap>::value_type spring_constant) |
| { |
| return kamada_kawai_spring_layout(g, position, weight, topology, edge_or_side_length, |
| done, spring_constant, |
| get(vertex_index, g)); |
| } |
| |
| /** |
| * \overload |
| */ |
| template<typename Topology, typename Graph, typename PositionMap, typename WeightMap, |
| typename T, bool EdgeOrSideLength, typename Done> |
| bool |
| kamada_kawai_spring_layout( |
| const Graph& g, |
| PositionMap position, |
| WeightMap weight, |
| const Topology& topology, |
| detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length, |
| Done done) |
| { |
| typedef typename property_traits<WeightMap>::value_type weight_type; |
| return kamada_kawai_spring_layout(g, position, weight, topology, edge_or_side_length, |
| done, weight_type(1)); |
| } |
| |
| /** |
| * \overload |
| */ |
| template<typename Topology, typename Graph, typename PositionMap, typename WeightMap, |
| typename T, bool EdgeOrSideLength> |
| bool |
| kamada_kawai_spring_layout( |
| const Graph& g, |
| PositionMap position, |
| WeightMap weight, |
| const Topology& topology, |
| detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length) |
| { |
| typedef typename property_traits<WeightMap>::value_type weight_type; |
| return kamada_kawai_spring_layout(g, position, weight, topology, edge_or_side_length, |
| layout_tolerance<weight_type>(), |
| weight_type(1.0), |
| get(vertex_index, g)); |
| } |
| } // end namespace boost |
| |
| #endif // BOOST_GRAPH_KAMADA_KAWAI_SPRING_LAYOUT_HPP |