williamr@2: // Copyright 2004 The Trustees of Indiana University. williamr@2: williamr@2: // Distributed under the Boost Software License, Version 1.0. williamr@2: // (See accompanying file LICENSE_1_0.txt or copy at williamr@2: // http://www.boost.org/LICENSE_1_0.txt) williamr@2: williamr@2: // Authors: Douglas Gregor williamr@2: // Andrew Lumsdaine williamr@2: #ifndef BOOST_GRAPH_KAMADA_KAWAI_SPRING_LAYOUT_HPP williamr@2: #define BOOST_GRAPH_KAMADA_KAWAI_SPRING_LAYOUT_HPP williamr@2: williamr@2: #include williamr@2: #include williamr@2: #include williamr@2: #include williamr@2: #include williamr@2: #include williamr@2: #include williamr@2: #include williamr@2: williamr@2: namespace boost { williamr@2: namespace detail { namespace graph { williamr@2: /** williamr@2: * Denotes an edge or display area side length used to scale a williamr@2: * Kamada-Kawai drawing. williamr@2: */ williamr@2: template williamr@2: struct edge_or_side williamr@2: { williamr@2: explicit edge_or_side(T value) : value(value) {} williamr@2: williamr@2: T value; williamr@2: }; williamr@2: williamr@2: /** williamr@2: * Compute the edge length from an edge length. This is trivial. williamr@2: */ williamr@2: template williamr@2: T compute_edge_length(const Graph&, DistanceMap, IndexMap, williamr@2: edge_or_side length) williamr@2: { return length.value; } williamr@2: williamr@2: /** williamr@2: * Compute the edge length based on the display area side williamr@2: length. We do this by dividing the side length by the largest williamr@2: shortest distance between any two vertices in the graph. williamr@2: */ williamr@2: template williamr@2: T williamr@2: compute_edge_length(const Graph& g, DistanceMap distance, IndexMap index, williamr@2: edge_or_side length) williamr@2: { williamr@2: T result(0); williamr@2: williamr@2: typedef typename graph_traits::vertex_iterator vertex_iterator; williamr@2: williamr@2: for (vertex_iterator ui = vertices(g).first, end = vertices(g).second; williamr@2: ui != end; ++ui) { williamr@2: vertex_iterator vi = ui; williamr@2: for (++vi; vi != end; ++vi) { williamr@2: T dij = distance[get(index, *ui)][get(index, *vi)]; williamr@2: if (dij > result) result = dij; williamr@2: } williamr@2: } williamr@2: return length.value / result; williamr@2: } williamr@2: williamr@2: /** williamr@2: * Implementation of the Kamada-Kawai spring layout algorithm. williamr@2: */ williamr@2: template williamr@2: struct kamada_kawai_spring_layout_impl williamr@2: { williamr@2: typedef typename property_traits::value_type weight_type; williamr@2: typedef std::pair deriv_type; williamr@2: typedef typename graph_traits::vertex_iterator vertex_iterator; williamr@2: typedef typename graph_traits::vertex_descriptor williamr@2: vertex_descriptor; williamr@2: williamr@2: kamada_kawai_spring_layout_impl( williamr@2: const Graph& g, williamr@2: PositionMap position, williamr@2: WeightMap weight, williamr@2: EdgeOrSideLength edge_or_side_length, williamr@2: Done done, williamr@2: weight_type spring_constant, williamr@2: VertexIndexMap index, williamr@2: DistanceMatrix distance, williamr@2: SpringStrengthMatrix spring_strength, williamr@2: PartialDerivativeMap partial_derivatives) williamr@2: : g(g), position(position), weight(weight), williamr@2: edge_or_side_length(edge_or_side_length), done(done), williamr@2: spring_constant(spring_constant), index(index), distance(distance), williamr@2: spring_strength(spring_strength), williamr@2: partial_derivatives(partial_derivatives) {} williamr@2: williamr@2: // Compute contribution of vertex i to the first partial williamr@2: // derivatives (dE/dx_m, dE/dy_m) (for vertex m) williamr@2: deriv_type williamr@2: compute_partial_derivative(vertex_descriptor m, vertex_descriptor i) williamr@2: { williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE williamr@2: using std::sqrt; williamr@2: #endif // BOOST_NO_STDC_NAMESPACE williamr@2: williamr@2: deriv_type result(0, 0); williamr@2: if (i != m) { williamr@2: weight_type x_diff = position[m].x - position[i].x; williamr@2: weight_type y_diff = position[m].y - position[i].y; williamr@2: weight_type dist = sqrt(x_diff * x_diff + y_diff * y_diff); williamr@2: result.first = spring_strength[get(index, m)][get(index, i)] williamr@2: * (x_diff - distance[get(index, m)][get(index, i)]*x_diff/dist); williamr@2: result.second = spring_strength[get(index, m)][get(index, i)] williamr@2: * (y_diff - distance[get(index, m)][get(index, i)]*y_diff/dist); williamr@2: } williamr@2: williamr@2: return result; williamr@2: } williamr@2: williamr@2: // Compute partial derivatives dE/dx_m and dE/dy_m williamr@2: deriv_type williamr@2: compute_partial_derivatives(vertex_descriptor m) williamr@2: { williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE williamr@2: using std::sqrt; williamr@2: #endif // BOOST_NO_STDC_NAMESPACE williamr@2: williamr@2: deriv_type result(0, 0); williamr@2: williamr@2: // TBD: looks like an accumulate to me williamr@2: std::pair verts = vertices(g); williamr@2: for (/* no init */; verts.first != verts.second; ++verts.first) { williamr@2: vertex_descriptor i = *verts.first; williamr@2: deriv_type deriv = compute_partial_derivative(m, i); williamr@2: result.first += deriv.first; williamr@2: result.second += deriv.second; williamr@2: } williamr@2: williamr@2: return result; williamr@2: } williamr@2: williamr@2: // The actual Kamada-Kawai spring layout algorithm implementation williamr@2: bool run() williamr@2: { williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE williamr@2: using std::sqrt; williamr@2: #endif // BOOST_NO_STDC_NAMESPACE williamr@2: williamr@2: // Compute d_{ij} and place it in the distance matrix williamr@2: if (!johnson_all_pairs_shortest_paths(g, distance, index, weight, williamr@2: weight_type(0))) williamr@2: return false; williamr@2: williamr@2: // Compute L based on side length (if needed), or retrieve L williamr@2: weight_type edge_length = williamr@2: detail::graph::compute_edge_length(g, distance, index, williamr@2: edge_or_side_length); williamr@2: williamr@2: // Compute l_{ij} and k_{ij} williamr@2: const weight_type K = spring_constant; williamr@2: vertex_iterator ui, end = vertices(g).second; williamr@2: for (ui = vertices(g).first; ui != end; ++ui) { williamr@2: vertex_iterator vi = ui; williamr@2: for (++vi; vi != end; ++vi) { williamr@2: weight_type dij = distance[get(index, *ui)][get(index, *vi)]; williamr@2: if (dij == (std::numeric_limits::max)()) williamr@2: return false; williamr@2: distance[get(index, *ui)][get(index, *vi)] = edge_length * dij; williamr@2: distance[get(index, *vi)][get(index, *ui)] = edge_length * dij; williamr@2: spring_strength[get(index, *ui)][get(index, *vi)] = K/(dij*dij); williamr@2: spring_strength[get(index, *vi)][get(index, *ui)] = K/(dij*dij); williamr@2: } williamr@2: } williamr@2: williamr@2: // Compute Delta_i and find max williamr@2: vertex_descriptor p = *vertices(g).first; williamr@2: weight_type delta_p(0); williamr@2: williamr@2: for (ui = vertices(g).first; ui != end; ++ui) { williamr@2: deriv_type deriv = compute_partial_derivatives(*ui); williamr@2: put(partial_derivatives, *ui, deriv); williamr@2: williamr@2: weight_type delta = williamr@2: sqrt(deriv.first*deriv.first + deriv.second*deriv.second); williamr@2: williamr@2: if (delta > delta_p) { williamr@2: p = *ui; williamr@2: delta_p = delta; williamr@2: } williamr@2: } williamr@2: williamr@2: while (!done(delta_p, p, g, true)) { williamr@2: // The contribution p makes to the partial derivatives of williamr@2: // each vertex. Computing this (at O(n) cost) allows us to williamr@2: // update the delta_i values in O(n) time instead of O(n^2) williamr@2: // time. williamr@2: std::vector p_partials(num_vertices(g)); williamr@2: for (ui = vertices(g).first; ui != end; ++ui) { williamr@2: vertex_descriptor i = *ui; williamr@2: p_partials[get(index, i)] = compute_partial_derivative(i, p); williamr@2: } williamr@2: williamr@2: do { williamr@2: // Compute the 4 elements of the Jacobian williamr@2: weight_type dE_dx_dx = 0, dE_dx_dy = 0, dE_dy_dx = 0, dE_dy_dy = 0; williamr@2: for (ui = vertices(g).first; ui != end; ++ui) { williamr@2: vertex_descriptor i = *ui; williamr@2: if (i != p) { williamr@2: weight_type x_diff = position[p].x - position[i].x; williamr@2: weight_type y_diff = position[p].y - position[i].y; williamr@2: weight_type dist = sqrt(x_diff * x_diff + y_diff * y_diff); williamr@2: weight_type dist_cubed = dist * dist * dist; williamr@2: weight_type k_mi = spring_strength[get(index,p)][get(index,i)]; williamr@2: weight_type l_mi = distance[get(index, p)][get(index, i)]; williamr@2: dE_dx_dx += k_mi * (1 - (l_mi * y_diff * y_diff)/dist_cubed); williamr@2: dE_dx_dy += k_mi * l_mi * x_diff * y_diff / dist_cubed; williamr@2: dE_dy_dx += k_mi * l_mi * x_diff * y_diff / dist_cubed; williamr@2: dE_dy_dy += k_mi * (1 - (l_mi * x_diff * x_diff)/dist_cubed); williamr@2: } williamr@2: } williamr@2: williamr@2: // Solve for delta_x and delta_y williamr@2: weight_type dE_dx = get(partial_derivatives, p).first; williamr@2: weight_type dE_dy = get(partial_derivatives, p).second; williamr@2: williamr@2: weight_type delta_x = williamr@2: (dE_dx_dy * dE_dy - dE_dy_dy * dE_dx) williamr@2: / (dE_dx_dx * dE_dy_dy - dE_dx_dy * dE_dy_dx); williamr@2: williamr@2: weight_type delta_y = williamr@2: (dE_dx_dx * dE_dy - dE_dy_dx * dE_dx) williamr@2: / (dE_dy_dx * dE_dx_dy - dE_dx_dx * dE_dy_dy); williamr@2: williamr@2: williamr@2: // Move p by (delta_x, delta_y) williamr@2: position[p].x += delta_x; williamr@2: position[p].y += delta_y; williamr@2: williamr@2: // Recompute partial derivatives and delta_p williamr@2: deriv_type deriv = compute_partial_derivatives(p); williamr@2: put(partial_derivatives, p, deriv); williamr@2: williamr@2: delta_p = williamr@2: sqrt(deriv.first*deriv.first + deriv.second*deriv.second); williamr@2: } while (!done(delta_p, p, g, false)); williamr@2: williamr@2: // Select new p by updating each partial derivative and delta williamr@2: vertex_descriptor old_p = p; williamr@2: for (ui = vertices(g).first; ui != end; ++ui) { williamr@2: deriv_type old_deriv_p = p_partials[get(index, *ui)]; williamr@2: deriv_type old_p_partial = williamr@2: compute_partial_derivative(*ui, old_p); williamr@2: deriv_type deriv = get(partial_derivatives, *ui); williamr@2: williamr@2: deriv.first += old_p_partial.first - old_deriv_p.first; williamr@2: deriv.second += old_p_partial.second - old_deriv_p.second; williamr@2: williamr@2: put(partial_derivatives, *ui, deriv); williamr@2: weight_type delta = williamr@2: sqrt(deriv.first*deriv.first + deriv.second*deriv.second); williamr@2: williamr@2: if (delta > delta_p) { williamr@2: p = *ui; williamr@2: delta_p = delta; williamr@2: } williamr@2: } williamr@2: } williamr@2: williamr@2: return true; williamr@2: } williamr@2: williamr@2: const Graph& g; williamr@2: PositionMap position; williamr@2: WeightMap weight; williamr@2: EdgeOrSideLength edge_or_side_length; williamr@2: Done done; williamr@2: weight_type spring_constant; williamr@2: VertexIndexMap index; williamr@2: DistanceMatrix distance; williamr@2: SpringStrengthMatrix spring_strength; williamr@2: PartialDerivativeMap partial_derivatives; williamr@2: }; williamr@2: } } // end namespace detail::graph williamr@2: williamr@2: /// States that the given quantity is an edge length. williamr@2: template williamr@2: inline detail::graph::edge_or_side williamr@2: edge_length(T x) williamr@2: { return detail::graph::edge_or_side(x); } williamr@2: williamr@2: /// States that the given quantity is a display area side length. williamr@2: template williamr@2: inline detail::graph::edge_or_side williamr@2: side_length(T x) williamr@2: { return detail::graph::edge_or_side(x); } williamr@2: williamr@2: /** williamr@2: * \brief Determines when to terminate layout of a particular graph based williamr@2: * on a given relative tolerance. williamr@2: */ williamr@2: template williamr@2: struct layout_tolerance williamr@2: { williamr@2: layout_tolerance(const T& tolerance = T(0.001)) williamr@2: : tolerance(tolerance), last_energy((std::numeric_limits::max)()), williamr@2: last_local_energy((std::numeric_limits::max)()) { } williamr@2: williamr@2: template williamr@2: bool williamr@2: operator()(T delta_p, williamr@2: typename boost::graph_traits::vertex_descriptor p, williamr@2: const Graph& g, williamr@2: bool global) williamr@2: { williamr@2: if (global) { williamr@2: if (last_energy == (std::numeric_limits::max)()) { williamr@2: last_energy = delta_p; williamr@2: return false; williamr@2: } williamr@2: williamr@2: T diff = last_energy - delta_p; williamr@2: if (diff < T(0)) diff = -diff; williamr@2: bool done = (delta_p == T(0) || diff / last_energy < tolerance); williamr@2: last_energy = delta_p; williamr@2: return done; williamr@2: } else { williamr@2: if (last_local_energy == (std::numeric_limits::max)()) { williamr@2: last_local_energy = delta_p; williamr@2: return delta_p == T(0); williamr@2: } williamr@2: williamr@2: T diff = last_local_energy - delta_p; williamr@2: bool done = (delta_p == T(0) || (diff / last_local_energy) < tolerance); williamr@2: last_local_energy = delta_p; williamr@2: return done; williamr@2: } williamr@2: } williamr@2: williamr@2: private: williamr@2: T tolerance; williamr@2: T last_energy; williamr@2: T last_local_energy; williamr@2: }; williamr@2: williamr@2: /** \brief Kamada-Kawai spring layout for undirected graphs. williamr@2: * williamr@2: * This algorithm performs graph layout (in two dimensions) for williamr@2: * connected, undirected graphs. It operates by relating the layout williamr@2: * of graphs to a dynamic spring system and minimizing the energy williamr@2: * within that system. The strength of a spring between two vertices williamr@2: * is inversely proportional to the square of the shortest distance williamr@2: * (in graph terms) between those two vertices. Essentially, williamr@2: * vertices that are closer in the graph-theoretic sense (i.e., by williamr@2: * following edges) will have stronger springs and will therefore be williamr@2: * placed closer together. williamr@2: * williamr@2: * Prior to invoking this algorithm, it is recommended that the williamr@2: * vertices be placed along the vertices of a regular n-sided williamr@2: * polygon. williamr@2: * williamr@2: * \param g (IN) must be a model of Vertex List Graph, Edge List williamr@2: * Graph, and Incidence Graph and must be undirected. williamr@2: * williamr@2: * \param position (OUT) must be a model of Lvalue Property Map, williamr@2: * where the value type is a class containing fields @c x and @c y williamr@2: * that will be set to the @c x and @c y coordinates of each vertex. williamr@2: * williamr@2: * \param weight (IN) must be a model of Readable Property Map, williamr@2: * which provides the weight of each edge in the graph @p g. williamr@2: * williamr@2: * \param edge_or_side_length (IN) provides either the unit length williamr@2: * @c e of an edge in the layout or the length of a side @c s of the williamr@2: * display area, and must be either @c boost::edge_length(e) or @c williamr@2: * boost::side_length(s), respectively. williamr@2: * williamr@2: * \param done (IN) is a 4-argument function object that is passed williamr@2: * the current value of delta_p (i.e., the energy of vertex @p p), williamr@2: * the vertex @p p, the graph @p g, and a boolean flag indicating williamr@2: * whether @p delta_p is the maximum energy in the system (when @c williamr@2: * true) or the energy of the vertex being moved. Defaults to @c williamr@2: * layout_tolerance instantiated over the value type of the weight williamr@2: * map. williamr@2: * williamr@2: * \param spring_constant (IN) is the constant multiplied by each williamr@2: * spring's strength. Larger values create systems with more energy williamr@2: * that can take longer to stabilize; smaller values create systems williamr@2: * with less energy that stabilize quickly but do not necessarily williamr@2: * result in pleasing layouts. The default value is 1. williamr@2: * williamr@2: * \param index (IN) is a mapping from vertices to index values williamr@2: * between 0 and @c num_vertices(g). The default is @c williamr@2: * get(vertex_index,g). williamr@2: * williamr@2: * \param distance (UTIL/OUT) will be used to store the distance williamr@2: * from every vertex to every other vertex, which is computed in the williamr@2: * first stages of the algorithm. This value's type must be a model williamr@2: * of BasicMatrix with value type equal to the value type of the williamr@2: * weight map. The default is a a vector of vectors. williamr@2: * williamr@2: * \param spring_strength (UTIL/OUT) will be used to store the williamr@2: * strength of the spring between every pair of vertices. This williamr@2: * value's type must be a model of BasicMatrix with value type equal williamr@2: * to the value type of the weight map. The default is a a vector of williamr@2: * vectors. williamr@2: * williamr@2: * \param partial_derivatives (UTIL) will be used to store the williamr@2: * partial derivates of each vertex with respect to the @c x and @c williamr@2: * y coordinates. This must be a Read/Write Property Map whose value williamr@2: * type is a pair with both types equivalent to the value type of williamr@2: * the weight map. The default is an iterator property map. williamr@2: * williamr@2: * \returns @c true if layout was successful or @c false if a williamr@2: * negative weight cycle was detected. williamr@2: */ williamr@2: template williamr@2: bool williamr@2: kamada_kawai_spring_layout( williamr@2: const Graph& g, williamr@2: PositionMap position, williamr@2: WeightMap weight, williamr@2: detail::graph::edge_or_side edge_or_side_length, williamr@2: Done done, williamr@2: typename property_traits::value_type spring_constant, williamr@2: VertexIndexMap index, williamr@2: DistanceMatrix distance, williamr@2: SpringStrengthMatrix spring_strength, williamr@2: PartialDerivativeMap partial_derivatives) williamr@2: { williamr@2: BOOST_STATIC_ASSERT((is_convertible< williamr@2: typename graph_traits::directed_category*, williamr@2: undirected_tag* williamr@2: >::value)); williamr@2: williamr@2: detail::graph::kamada_kawai_spring_layout_impl< williamr@2: Graph, PositionMap, WeightMap, williamr@2: detail::graph::edge_or_side, Done, VertexIndexMap, williamr@2: DistanceMatrix, SpringStrengthMatrix, PartialDerivativeMap> williamr@2: alg(g, position, weight, edge_or_side_length, done, spring_constant, williamr@2: index, distance, spring_strength, partial_derivatives); williamr@2: return alg.run(); williamr@2: } williamr@2: williamr@2: /** williamr@2: * \overload williamr@2: */ williamr@2: template williamr@2: bool williamr@2: kamada_kawai_spring_layout( williamr@2: const Graph& g, williamr@2: PositionMap position, williamr@2: WeightMap weight, williamr@2: detail::graph::edge_or_side edge_or_side_length, williamr@2: Done done, williamr@2: typename property_traits::value_type spring_constant, williamr@2: VertexIndexMap index) williamr@2: { williamr@2: typedef typename property_traits::value_type weight_type; williamr@2: williamr@2: typename graph_traits::vertices_size_type n = num_vertices(g); williamr@2: typedef std::vector weight_vec; williamr@2: williamr@2: std::vector distance(n, weight_vec(n)); williamr@2: std::vector spring_strength(n, weight_vec(n)); williamr@2: std::vector > partial_derivatives(n); williamr@2: williamr@2: return williamr@2: kamada_kawai_spring_layout( williamr@2: g, position, weight, edge_or_side_length, done, spring_constant, index, williamr@2: distance.begin(), williamr@2: spring_strength.begin(), williamr@2: make_iterator_property_map(partial_derivatives.begin(), index, williamr@2: std::pair())); williamr@2: } williamr@2: williamr@2: /** williamr@2: * \overload williamr@2: */ williamr@2: template williamr@2: bool williamr@2: kamada_kawai_spring_layout( williamr@2: const Graph& g, williamr@2: PositionMap position, williamr@2: WeightMap weight, williamr@2: detail::graph::edge_or_side edge_or_side_length, williamr@2: Done done, williamr@2: typename property_traits::value_type spring_constant) williamr@2: { williamr@2: return kamada_kawai_spring_layout(g, position, weight, edge_or_side_length, williamr@2: done, spring_constant, williamr@2: get(vertex_index, g)); williamr@2: } williamr@2: williamr@2: /** williamr@2: * \overload williamr@2: */ williamr@2: template williamr@2: bool williamr@2: kamada_kawai_spring_layout( williamr@2: const Graph& g, williamr@2: PositionMap position, williamr@2: WeightMap weight, williamr@2: detail::graph::edge_or_side edge_or_side_length, williamr@2: Done done) williamr@2: { williamr@2: typedef typename property_traits::value_type weight_type; williamr@2: return kamada_kawai_spring_layout(g, position, weight, edge_or_side_length, williamr@2: done, weight_type(1)); williamr@2: } williamr@2: williamr@2: /** williamr@2: * \overload williamr@2: */ williamr@2: template williamr@2: bool williamr@2: kamada_kawai_spring_layout( williamr@2: const Graph& g, williamr@2: PositionMap position, williamr@2: WeightMap weight, williamr@2: detail::graph::edge_or_side edge_or_side_length) williamr@2: { williamr@2: typedef typename property_traits::value_type weight_type; williamr@2: return kamada_kawai_spring_layout(g, position, weight, edge_or_side_length, williamr@2: layout_tolerance(), williamr@2: weight_type(1.0), williamr@2: get(vertex_index, g)); williamr@2: } williamr@2: } // end namespace boost williamr@2: williamr@2: #endif // BOOST_GRAPH_KAMADA_KAWAI_SPRING_LAYOUT_HPP