1 // Copyright 2004 The Trustees of Indiana University.
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)
7 // Authors: Douglas Gregor
9 #ifndef BOOST_GRAPH_KAMADA_KAWAI_SPRING_LAYOUT_HPP
10 #define BOOST_GRAPH_KAMADA_KAWAI_SPRING_LAYOUT_HPP
12 #include <boost/graph/graph_traits.hpp>
13 #include <boost/graph/topology.hpp>
14 #include <boost/graph/iteration_macros.hpp>
15 #include <boost/graph/johnson_all_pairs_shortest.hpp>
16 #include <boost/type_traits/is_convertible.hpp>
21 #include <boost/limits.hpp>
22 #include <boost/config/no_tr1/cmath.hpp>
25 namespace detail { namespace graph {
27 * Denotes an edge or display area side length used to scale a
28 * Kamada-Kawai drawing.
30 template<bool Edge, typename T>
33 explicit edge_or_side(T value) : value(value) {}
39 * Compute the edge length from an edge length. This is trivial.
41 template<typename Graph, typename DistanceMap, typename IndexMap,
43 T compute_edge_length(const Graph&, DistanceMap, IndexMap,
44 edge_or_side<true, T> length)
45 { return length.value; }
48 * Compute the edge length based on the display area side
49 length. We do this by dividing the side length by the largest
50 shortest distance between any two vertices in the graph.
52 template<typename Graph, typename DistanceMap, typename IndexMap,
55 compute_edge_length(const Graph& g, DistanceMap distance, IndexMap index,
56 edge_or_side<false, T> length)
60 typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
62 for (vertex_iterator ui = vertices(g).first, end = vertices(g).second;
64 vertex_iterator vi = ui;
65 for (++vi; vi != end; ++vi) {
66 T dij = distance[get(index, *ui)][get(index, *vi)];
67 if (dij > result) result = dij;
70 return length.value / result;
74 * Dense linear solver for fixed-size matrices.
76 template <std::size_t Size>
77 struct linear_solver {
78 // Indices in mat are (row, column)
79 // template <typename Vec>
80 // static Vec solve(double mat[Size][Size], Vec rhs);
84 struct linear_solver<1> {
85 template <typename Vec>
86 static Vec solve(double mat[1][1], Vec rhs) {
87 return rhs / mat[0][0];
91 // These are from http://en.wikipedia.org/wiki/Cramer%27s_rule
93 struct linear_solver<2> {
94 template <typename Vec>
95 static Vec solve(double mat[2][2], Vec rhs) {
96 double denom = mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1];
97 double x_num = rhs[0] * mat[1][1] - rhs[1] * mat[0][1];
98 double y_num = mat[0][0] * rhs[1] - mat[1][0] * rhs[0] ;
100 result[0] = x_num / denom;
101 result[1] = y_num / denom;
107 struct linear_solver<3> {
108 template <typename Vec>
109 static Vec solve(double mat[3][3], Vec rhs) {
110 double denom = mat[0][0] * (mat[1][1] * mat[2][2] - mat[2][1] * mat[1][2])
111 - mat[1][0] * (mat[0][1] * mat[2][2] - mat[2][1] * mat[0][2])
112 + mat[2][0] * (mat[0][1] * mat[1][2] - mat[1][1] * mat[0][2]);
113 double x_num = rhs[0] * (mat[1][1] * mat[2][2] - mat[2][1] * mat[1][2])
114 - rhs[1] * (mat[0][1] * mat[2][2] - mat[2][1] * mat[0][2])
115 + rhs[2] * (mat[0][1] * mat[1][2] - mat[1][1] * mat[0][2]);
116 double y_num = mat[0][0] * (rhs[1] * mat[2][2] - rhs[2] * mat[1][2])
117 - mat[1][0] * (rhs[0] * mat[2][2] - rhs[2] * mat[0][2])
118 + mat[2][0] * (rhs[0] * mat[1][2] - rhs[1] * mat[0][2]);
119 double z_num = mat[0][0] * (mat[1][1] * rhs[2] - mat[2][1] * rhs[1] )
120 - mat[1][0] * (mat[0][1] * rhs[2] - mat[2][1] * rhs[0] )
121 + mat[2][0] * (mat[0][1] * rhs[1] - mat[1][1] * rhs[0] );
123 result[0] = x_num / denom;
124 result[1] = y_num / denom;
125 result[2] = z_num / denom;
131 * Implementation of the Kamada-Kawai spring layout algorithm.
133 template<typename Topology, typename Graph, typename PositionMap, typename WeightMap,
134 typename EdgeOrSideLength, typename Done,
135 typename VertexIndexMap, typename DistanceMatrix,
136 typename SpringStrengthMatrix, typename PartialDerivativeMap>
137 struct kamada_kawai_spring_layout_impl
139 typedef typename property_traits<WeightMap>::value_type weight_type;
140 typedef typename Topology::point_type Point;
141 typedef typename Topology::point_difference_type point_difference_type;
142 typedef point_difference_type deriv_type;
143 typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
144 typedef typename graph_traits<Graph>::vertex_descriptor
147 kamada_kawai_spring_layout_impl(
148 const Topology& topology,
150 PositionMap position,
152 EdgeOrSideLength edge_or_side_length,
154 weight_type spring_constant,
155 VertexIndexMap index,
156 DistanceMatrix distance,
157 SpringStrengthMatrix spring_strength,
158 PartialDerivativeMap partial_derivatives)
159 : topology(topology), g(g), position(position), weight(weight),
160 edge_or_side_length(edge_or_side_length), done(done),
161 spring_constant(spring_constant), index(index), distance(distance),
162 spring_strength(spring_strength),
163 partial_derivatives(partial_derivatives) {}
165 // Compute contribution of vertex i to the first partial
166 // derivatives (dE/dx_m, dE/dy_m) (for vertex m)
168 compute_partial_derivative(vertex_descriptor m, vertex_descriptor i)
170 #ifndef BOOST_NO_STDC_NAMESPACE
172 #endif // BOOST_NO_STDC_NAMESPACE
176 point_difference_type diff = topology.difference(position[m], position[i]);
177 weight_type dist = topology.norm(diff);
178 result = spring_strength[get(index, m)][get(index, i)]
179 * (diff - distance[get(index, m)][get(index, i)]/dist*diff);
185 // Compute partial derivatives dE/dx_m and dE/dy_m
187 compute_partial_derivatives(vertex_descriptor m)
189 #ifndef BOOST_NO_STDC_NAMESPACE
191 #endif // BOOST_NO_STDC_NAMESPACE
195 // TBD: looks like an accumulate to me
196 BGL_FORALL_VERTICES_T(i, g, Graph) {
197 deriv_type deriv = compute_partial_derivative(m, i);
204 // The actual Kamada-Kawai spring layout algorithm implementation
207 #ifndef BOOST_NO_STDC_NAMESPACE
209 #endif // BOOST_NO_STDC_NAMESPACE
211 // Compute d_{ij} and place it in the distance matrix
212 if (!johnson_all_pairs_shortest_paths(g, distance, index, weight,
216 // Compute L based on side length (if needed), or retrieve L
217 weight_type edge_length =
218 detail::graph::compute_edge_length(g, distance, index,
219 edge_or_side_length);
221 // std::cerr << "edge_length = " << edge_length << std::endl;
223 // Compute l_{ij} and k_{ij}
224 const weight_type K = spring_constant;
225 vertex_iterator ui, end;
226 for (ui = vertices(g).first, end = vertices(g).second; ui != end; ++ui) {
227 vertex_iterator vi = ui;
228 for (++vi; vi != end; ++vi) {
229 weight_type dij = distance[get(index, *ui)][get(index, *vi)];
230 if (dij == (std::numeric_limits<weight_type>::max)())
232 distance[get(index, *ui)][get(index, *vi)] = edge_length * dij;
233 distance[get(index, *vi)][get(index, *ui)] = edge_length * dij;
234 spring_strength[get(index, *ui)][get(index, *vi)] = K/(dij*dij);
235 spring_strength[get(index, *vi)][get(index, *ui)] = K/(dij*dij);
239 // Compute Delta_i and find max
240 vertex_descriptor p = *vertices(g).first;
241 weight_type delta_p(0);
243 for (ui = vertices(g).first, end = vertices(g).second; ui != end; ++ui) {
244 deriv_type deriv = compute_partial_derivatives(*ui);
245 put(partial_derivatives, *ui, deriv);
247 weight_type delta = topology.norm(deriv);
249 if (delta > delta_p) {
255 while (!done(delta_p, p, g, true)) {
256 // The contribution p makes to the partial derivatives of
257 // each vertex. Computing this (at O(n) cost) allows us to
258 // update the delta_i values in O(n) time instead of O(n^2)
260 std::vector<deriv_type> p_partials(num_vertices(g));
261 for (ui = vertices(g).first, end = vertices(g).second; ui != end; ++ui) {
262 vertex_descriptor i = *ui;
263 p_partials[get(index, i)] = compute_partial_derivative(i, p);
267 // For debugging, compute the energy value E
269 for (ui = vertices(g).first, end = vertices(g).second; ui != end; ++ui) {
270 vertex_iterator vi = ui;
271 for (++vi; vi != end; ++vi) {
272 double dist = topology.distance(position[*ui], position[*vi]);
273 weight_type k_ij = spring_strength[get(index,*ui)][get(index,*vi)];
274 weight_type l_ij = distance[get(index, *ui)][get(index, *vi)];
275 E += .5 * k_ij * (dist - l_ij) * (dist - l_ij);
278 // std::cerr << "E = " << E << std::endl;
280 // Compute the elements of the Jacobian
282 // http://www.cs.panam.edu/~rfowler/papers/1994_kumar_fowler_A_Spring_UTPACSTR.pdf
283 // with the bugs fixed in the off-diagonal case
284 weight_type dE_d_d[Point::dimensions][Point::dimensions];
285 for (std::size_t i = 0; i < Point::dimensions; ++i)
286 for (std::size_t j = 0; j < Point::dimensions; ++j)
288 for (ui = vertices(g).first, end = vertices(g).second; ui != end; ++ui) {
289 vertex_descriptor i = *ui;
291 point_difference_type diff = topology.difference(position[p], position[i]);
292 weight_type dist = topology.norm(diff);
293 weight_type dist_squared = dist * dist;
294 weight_type inv_dist_cubed = 1. / (dist_squared * dist);
295 weight_type k_mi = spring_strength[get(index,p)][get(index,i)];
296 weight_type l_mi = distance[get(index, p)][get(index, i)];
297 for (std::size_t i = 0; i < Point::dimensions; ++i) {
298 for (std::size_t j = 0; j < Point::dimensions; ++j) {
300 dE_d_d[i][i] += k_mi * (1 + (l_mi * (diff[i] * diff[i] - dist_squared) * inv_dist_cubed));
302 dE_d_d[i][j] += k_mi * l_mi * diff[i] * diff[j] * inv_dist_cubed;
303 // dE_d_d[i][j] += k_mi * l_mi * sqrt(hypot(diff[i], diff[j])) * inv_dist_cubed;
310 deriv_type dE_d = get(partial_derivatives, p);
312 // Solve dE_d_d * delta = -dE_d to get delta
313 point_difference_type delta = -linear_solver<Point::dimensions>::solve(dE_d_d, dE_d);
316 position[p] = topology.adjust(position[p], delta);
318 // Recompute partial derivatives and delta_p
319 deriv_type deriv = compute_partial_derivatives(p);
320 put(partial_derivatives, p, deriv);
322 delta_p = topology.norm(deriv);
323 } while (!done(delta_p, p, g, false));
325 // Select new p by updating each partial derivative and delta
326 vertex_descriptor old_p = p;
327 for (ui = vertices(g).first, end = vertices(g).second; ui != end; ++ui) {
328 deriv_type old_deriv_p = p_partials[get(index, *ui)];
329 deriv_type old_p_partial =
330 compute_partial_derivative(*ui, old_p);
331 deriv_type deriv = get(partial_derivatives, *ui);
333 deriv += old_p_partial - old_deriv_p;
335 put(partial_derivatives, *ui, deriv);
336 weight_type delta = topology.norm(deriv);
338 if (delta > delta_p) {
348 const Topology& topology;
350 PositionMap position;
352 EdgeOrSideLength edge_or_side_length;
354 weight_type spring_constant;
355 VertexIndexMap index;
356 DistanceMatrix distance;
357 SpringStrengthMatrix spring_strength;
358 PartialDerivativeMap partial_derivatives;
360 } } // end namespace detail::graph
362 /// States that the given quantity is an edge length.
364 inline detail::graph::edge_or_side<true, T>
366 { return detail::graph::edge_or_side<true, T>(x); }
368 /// States that the given quantity is a display area side length.
370 inline detail::graph::edge_or_side<false, T>
372 { return detail::graph::edge_or_side<false, T>(x); }
375 * \brief Determines when to terminate layout of a particular graph based
376 * on a given relative tolerance.
378 template<typename T = double>
379 struct layout_tolerance
381 layout_tolerance(const T& tolerance = T(0.001))
382 : tolerance(tolerance), last_energy((std::numeric_limits<T>::max)()),
383 last_local_energy((std::numeric_limits<T>::max)()) { }
385 template<typename Graph>
387 operator()(T delta_p,
388 typename boost::graph_traits<Graph>::vertex_descriptor p,
393 if (last_energy == (std::numeric_limits<T>::max)()) {
394 last_energy = delta_p;
398 T diff = last_energy - delta_p;
399 if (diff < T(0)) diff = -diff;
400 bool done = (delta_p == T(0) || diff / last_energy < tolerance);
401 last_energy = delta_p;
404 if (last_local_energy == (std::numeric_limits<T>::max)()) {
405 last_local_energy = delta_p;
406 return delta_p == T(0);
409 T diff = last_local_energy - delta_p;
410 bool done = (delta_p == T(0) || (diff / last_local_energy) < tolerance);
411 last_local_energy = delta_p;
422 /** \brief Kamada-Kawai spring layout for undirected graphs.
424 * This algorithm performs graph layout (in two dimensions) for
425 * connected, undirected graphs. It operates by relating the layout
426 * of graphs to a dynamic spring system and minimizing the energy
427 * within that system. The strength of a spring between two vertices
428 * is inversely proportional to the square of the shortest distance
429 * (in graph terms) between those two vertices. Essentially,
430 * vertices that are closer in the graph-theoretic sense (i.e., by
431 * following edges) will have stronger springs and will therefore be
432 * placed closer together.
434 * Prior to invoking this algorithm, it is recommended that the
435 * vertices be placed along the vertices of a regular n-sided
438 * \param g (IN) must be a model of Vertex List Graph, Edge List
439 * Graph, and Incidence Graph and must be undirected.
441 * \param position (OUT) must be a model of Lvalue Property Map,
442 * where the value type is a class containing fields @c x and @c y
443 * that will be set to the @c x and @c y coordinates of each vertex.
445 * \param weight (IN) must be a model of Readable Property Map,
446 * which provides the weight of each edge in the graph @p g.
448 * \param topology (IN) must be a topology object (see topology.hpp),
449 * which provides operations on points and differences between them.
451 * \param edge_or_side_length (IN) provides either the unit length
452 * @c e of an edge in the layout or the length of a side @c s of the
453 * display area, and must be either @c boost::edge_length(e) or @c
454 * boost::side_length(s), respectively.
456 * \param done (IN) is a 4-argument function object that is passed
457 * the current value of delta_p (i.e., the energy of vertex @p p),
458 * the vertex @p p, the graph @p g, and a boolean flag indicating
459 * whether @p delta_p is the maximum energy in the system (when @c
460 * true) or the energy of the vertex being moved. Defaults to @c
461 * layout_tolerance instantiated over the value type of the weight
464 * \param spring_constant (IN) is the constant multiplied by each
465 * spring's strength. Larger values create systems with more energy
466 * that can take longer to stabilize; smaller values create systems
467 * with less energy that stabilize quickly but do not necessarily
468 * result in pleasing layouts. The default value is 1.
470 * \param index (IN) is a mapping from vertices to index values
471 * between 0 and @c num_vertices(g). The default is @c
472 * get(vertex_index,g).
474 * \param distance (UTIL/OUT) will be used to store the distance
475 * from every vertex to every other vertex, which is computed in the
476 * first stages of the algorithm. This value's type must be a model
477 * of BasicMatrix with value type equal to the value type of the
478 * weight map. The default is a vector of vectors.
480 * \param spring_strength (UTIL/OUT) will be used to store the
481 * strength of the spring between every pair of vertices. This
482 * value's type must be a model of BasicMatrix with value type equal
483 * to the value type of the weight map. The default is a vector of
486 * \param partial_derivatives (UTIL) will be used to store the
487 * partial derivates of each vertex with respect to the @c x and @c
488 * y coordinates. This must be a Read/Write Property Map whose value
489 * type is a pair with both types equivalent to the value type of
490 * the weight map. The default is an iterator property map.
492 * \returns @c true if layout was successful or @c false if a
493 * negative weight cycle was detected.
495 template<typename Topology, typename Graph, typename PositionMap, typename WeightMap,
496 typename T, bool EdgeOrSideLength, typename Done,
497 typename VertexIndexMap, typename DistanceMatrix,
498 typename SpringStrengthMatrix, typename PartialDerivativeMap>
500 kamada_kawai_spring_layout(
502 PositionMap position,
504 const Topology& topology,
505 detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length,
507 typename property_traits<WeightMap>::value_type spring_constant,
508 VertexIndexMap index,
509 DistanceMatrix distance,
510 SpringStrengthMatrix spring_strength,
511 PartialDerivativeMap partial_derivatives)
513 BOOST_STATIC_ASSERT((is_convertible<
514 typename graph_traits<Graph>::directed_category*,
518 detail::graph::kamada_kawai_spring_layout_impl<
519 Topology, Graph, PositionMap, WeightMap,
520 detail::graph::edge_or_side<EdgeOrSideLength, T>, Done, VertexIndexMap,
521 DistanceMatrix, SpringStrengthMatrix, PartialDerivativeMap>
522 alg(topology, g, position, weight, edge_or_side_length, done, spring_constant,
523 index, distance, spring_strength, partial_derivatives);
530 template<typename Topology, typename Graph, typename PositionMap, typename WeightMap,
531 typename T, bool EdgeOrSideLength, typename Done,
532 typename VertexIndexMap>
534 kamada_kawai_spring_layout(
536 PositionMap position,
538 const Topology& topology,
539 detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length,
541 typename property_traits<WeightMap>::value_type spring_constant,
542 VertexIndexMap index)
544 typedef typename property_traits<WeightMap>::value_type weight_type;
546 typename graph_traits<Graph>::vertices_size_type n = num_vertices(g);
547 typedef std::vector<weight_type> weight_vec;
549 std::vector<weight_vec> distance(n, weight_vec(n));
550 std::vector<weight_vec> spring_strength(n, weight_vec(n));
551 std::vector<typename Topology::point_difference_type> partial_derivatives(n);
554 kamada_kawai_spring_layout(
555 g, position, weight, topology, edge_or_side_length, done, spring_constant, index,
557 spring_strength.begin(),
558 make_iterator_property_map(partial_derivatives.begin(), index,
559 typename Topology::point_difference_type()));
565 template<typename Topology, typename Graph, typename PositionMap, typename WeightMap,
566 typename T, bool EdgeOrSideLength, typename Done>
568 kamada_kawai_spring_layout(
570 PositionMap position,
572 const Topology& topology,
573 detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length,
575 typename property_traits<WeightMap>::value_type spring_constant)
577 return kamada_kawai_spring_layout(g, position, weight, topology, edge_or_side_length,
578 done, spring_constant,
579 get(vertex_index, g));
585 template<typename Topology, typename Graph, typename PositionMap, typename WeightMap,
586 typename T, bool EdgeOrSideLength, typename Done>
588 kamada_kawai_spring_layout(
590 PositionMap position,
592 const Topology& topology,
593 detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length,
596 typedef typename property_traits<WeightMap>::value_type weight_type;
597 return kamada_kawai_spring_layout(g, position, weight, topology, edge_or_side_length,
598 done, weight_type(1));
604 template<typename Topology, typename Graph, typename PositionMap, typename WeightMap,
605 typename T, bool EdgeOrSideLength>
607 kamada_kawai_spring_layout(
609 PositionMap position,
611 const Topology& topology,
612 detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length)
614 typedef typename property_traits<WeightMap>::value_type weight_type;
615 return kamada_kawai_spring_layout(g, position, weight, topology, edge_or_side_length,
616 layout_tolerance<weight_type>(),
618 get(vertex_index, g));
620 } // end namespace boost
622 #endif // BOOST_GRAPH_KAMADA_KAWAI_SPRING_LAYOUT_HPP