1// Copyright 2004 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_KAMADA_KAWAI_SPRING_LAYOUT_HPP
10#define BOOST_GRAPH_KAMADA_KAWAI_SPRING_LAYOUT_HPP
11
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>
17#include <utility>
18#include <iterator>
19#include <vector>
20#include <iostream>
21#include <boost/limits.hpp>
22#include <boost/config/no_tr1/cmath.hpp>
23
24namespace boost {
25 namespace detail { namespace graph {
26 /**
27 * Denotes an edge or display area side length used to scale a
28 * Kamada-Kawai drawing.
29 */
30 template<bool Edge, typename T>
31 struct edge_or_side
32 {
33 explicit edge_or_side(T value) : value(value) {}
34
35 T value;
36 };
37
38 /**
39 * Compute the edge length from an edge length. This is trivial.
40 */
41 template<typename Graph, typename DistanceMap, typename IndexMap,
42 typename T>
43 T compute_edge_length(const Graph&, DistanceMap, IndexMap,
44 edge_or_side<true, T> length)
45 { return length.value; }
46
47 /**
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.
51 */
52 template<typename Graph, typename DistanceMap, typename IndexMap,
53 typename T>
54 T
55 compute_edge_length(const Graph& g, DistanceMap distance, IndexMap index,
56 edge_or_side<false, T> length)
57 {
58 T result(0);
59
60 typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
61
62 for (vertex_iterator ui = vertices(g).first, end = vertices(g).second;
63 ui != end; ++ui) {
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;
68 }
69 }
70 return length.value / result;
71 }
72
73 /**
74 * Dense linear solver for fixed-size matrices.
75 */
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);
81 };
82
83 template <>
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];
88 }
89 };
90
91 // These are from http://en.wikipedia.org/wiki/Cramer%27s_rule
92 template <>
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] ;
99 Vec result;
100 result[0] = x_num / denom;
101 result[1] = y_num / denom;
102 return result;
103 }
104 };
105
106 template <>
107 struct linear_solver<3> {
108 template <typename Vec>
109 static Vec solve(double mat[2][2], 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] );
122 Vec result;
123 result[0] = x_num / denom;
124 result[1] = y_num / denom;
125 result[2] = z_num / denom;
126 return result;
127 }
128 };
129
130 /**
131 * Implementation of the Kamada-Kawai spring layout algorithm.
132 */
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
138 {
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
145 vertex_descriptor;
146
147 kamada_kawai_spring_layout_impl(
148 const Topology& topology,
149 const Graph& g,
150 PositionMap position,
151 WeightMap weight,
152 EdgeOrSideLength edge_or_side_length,
153 Done done,
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) {}
164
165 // Compute contribution of vertex i to the first partial
166 // derivatives (dE/dx_m, dE/dy_m) (for vertex m)
167 deriv_type
168 compute_partial_derivative(vertex_descriptor m, vertex_descriptor i)
169 {
170#ifndef BOOST_NO_STDC_NAMESPACE
171 using std::sqrt;
172#endif // BOOST_NO_STDC_NAMESPACE
173
174 deriv_type result;
175 if (i != m) {
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);
180 }
181
182 return result;
183 }
184
185 // Compute partial derivatives dE/dx_m and dE/dy_m
186 deriv_type
187 compute_partial_derivatives(vertex_descriptor m)
188 {
189#ifndef BOOST_NO_STDC_NAMESPACE
190 using std::sqrt;
191#endif // BOOST_NO_STDC_NAMESPACE
192
193 deriv_type result;
194
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);
198 result += deriv;
199 }
200
201 return result;
202 }
203
204 // The actual Kamada-Kawai spring layout algorithm implementation
205 bool run()
206 {
207#ifndef BOOST_NO_STDC_NAMESPACE
208 using std::sqrt;
209#endif // BOOST_NO_STDC_NAMESPACE
210
211 // Compute d_{ij} and place it in the distance matrix
212 if (!johnson_all_pairs_shortest_paths(g, distance, index, weight,
213 weight_type(0)))
214 return false;
215
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);
220
221 // std::cerr << "edge_length = " << edge_length << std::endl;
222
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)())
231 return false;
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);
236 }
237 }
238
239 // Compute Delta_i and find max
240 vertex_descriptor p = *vertices(g).first;
241 weight_type delta_p(0);
242
243 for (ui = vertices(g).first, end = vertices(g).second; ui != end; ++ui) {
244 deriv_type deriv = compute_partial_derivatives(m: *ui);
245 put(partial_derivatives, *ui, deriv);
246
247 weight_type delta = topology.norm(deriv);
248
249 if (delta > delta_p) {
250 p = *ui;
251 delta_p = delta;
252 }
253 }
254
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)
259 // time.
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(m: i, i: p);
264 }
265
266 do {
267 // For debugging, compute the energy value E
268 double E = 0.;
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);
276 }
277 }
278 // std::cerr << "E = " << E << std::endl;
279
280 // Compute the elements of the Jacobian
281 // From
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)
287 dE_d_d[i][j] = 0.;
288 for (ui = vertices(g).first, end = vertices(g).second; ui != end; ++ui) {
289 vertex_descriptor i = *ui;
290 if (i != p) {
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) {
299 if (i == j) {
300 dE_d_d[i][i] += k_mi * (1 + (l_mi * (diff[i] * diff[i] - dist_squared) * inv_dist_cubed));
301 } else {
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;
304 }
305 }
306 }
307 }
308 }
309
310 deriv_type dE_d = get(partial_derivatives, p);
311
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);
314
315 // Move p by delta
316 position[p] = topology.adjust(position[p], delta);
317
318 // Recompute partial derivatives and delta_p
319 deriv_type deriv = compute_partial_derivatives(m: p);
320 put(partial_derivatives, p, deriv);
321
322 delta_p = topology.norm(deriv);
323 } while (!done(delta_p, p, g, false));
324
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(m: *ui, i: old_p);
331 deriv_type deriv = get(partial_derivatives, *ui);
332
333 deriv += old_p_partial - old_deriv_p;
334
335 put(partial_derivatives, *ui, deriv);
336 weight_type delta = topology.norm(deriv);
337
338 if (delta > delta_p) {
339 p = *ui;
340 delta_p = delta;
341 }
342 }
343 }
344
345 return true;
346 }
347
348 const Topology& topology;
349 const Graph& g;
350 PositionMap position;
351 WeightMap weight;
352 EdgeOrSideLength edge_or_side_length;
353 Done done;
354 weight_type spring_constant;
355 VertexIndexMap index;
356 DistanceMatrix distance;
357 SpringStrengthMatrix spring_strength;
358 PartialDerivativeMap partial_derivatives;
359 };
360 } } // end namespace detail::graph
361
362 /// States that the given quantity is an edge length.
363 template<typename T>
364 inline detail::graph::edge_or_side<true, T>
365 edge_length(T x)
366 { return detail::graph::edge_or_side<true, T>(x); }
367
368 /// States that the given quantity is a display area side length.
369 template<typename T>
370 inline detail::graph::edge_or_side<false, T>
371 side_length(T x)
372 { return detail::graph::edge_or_side<false, T>(x); }
373
374 /**
375 * \brief Determines when to terminate layout of a particular graph based
376 * on a given relative tolerance.
377 */
378 template<typename T = double>
379 struct layout_tolerance
380 {
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)()) { }
384
385 template<typename Graph>
386 bool
387 operator()(T delta_p,
388 typename boost::graph_traits<Graph>::vertex_descriptor p,
389 const Graph& g,
390 bool global)
391 {
392 if (global) {
393 if (last_energy == (std::numeric_limits<T>::max)()) {
394 last_energy = delta_p;
395 return false;
396 }
397
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;
402 return done;
403 } else {
404 if (last_local_energy == (std::numeric_limits<T>::max)()) {
405 last_local_energy = delta_p;
406 return delta_p == T(0);
407 }
408
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;
412 return done;
413 }
414 }
415
416 private:
417 T tolerance;
418 T last_energy;
419 T last_local_energy;
420 };
421
422 /** \brief Kamada-Kawai spring layout for undirected graphs.
423 *
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.
433 *
434 * Prior to invoking this algorithm, it is recommended that the
435 * vertices be placed along the vertices of a regular n-sided
436 * polygon.
437 *
438 * \param g (IN) must be a model of Vertex List Graph, Edge List
439 * Graph, and Incidence Graph and must be undirected.
440 *
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.
444 *
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.
447 *
448 * \param topology (IN) must be a topology object (see topology.hpp),
449 * which provides operations on points and differences between them.
450 *
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.
455 *
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
462 * map.
463 *
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.
469 *
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).
473 *
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.
479 *
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
484 * vectors.
485 *
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.
491 *
492 * \returns @c true if layout was successful or @c false if a
493 * negative weight cycle was detected.
494 */
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>
499 bool
500 kamada_kawai_spring_layout(
501 const Graph& g,
502 PositionMap position,
503 WeightMap weight,
504 const Topology& topology,
505 detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length,
506 Done done,
507 typename property_traits<WeightMap>::value_type spring_constant,
508 VertexIndexMap index,
509 DistanceMatrix distance,
510 SpringStrengthMatrix spring_strength,
511 PartialDerivativeMap partial_derivatives)
512 {
513 BOOST_STATIC_ASSERT((is_convertible<
514 typename graph_traits<Graph>::directed_category*,
515 undirected_tag*
516 >::value));
517
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);
524 return alg.run();
525 }
526
527 /**
528 * \overload
529 */
530 template<typename Topology, typename Graph, typename PositionMap, typename WeightMap,
531 typename T, bool EdgeOrSideLength, typename Done,
532 typename VertexIndexMap>
533 bool
534 kamada_kawai_spring_layout(
535 const Graph& g,
536 PositionMap position,
537 WeightMap weight,
538 const Topology& topology,
539 detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length,
540 Done done,
541 typename property_traits<WeightMap>::value_type spring_constant,
542 VertexIndexMap index)
543 {
544 typedef typename property_traits<WeightMap>::value_type weight_type;
545
546 typename graph_traits<Graph>::vertices_size_type n = num_vertices(g);
547 typedef std::vector<weight_type> weight_vec;
548
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);
552
553 return
554 kamada_kawai_spring_layout(
555 g, position, weight, topology, edge_or_side_length, done, spring_constant, index,
556 distance.begin(),
557 spring_strength.begin(),
558 make_iterator_property_map(partial_derivatives.begin(), index,
559 typename Topology::point_difference_type()));
560 }
561
562 /**
563 * \overload
564 */
565 template<typename Topology, typename Graph, typename PositionMap, typename WeightMap,
566 typename T, bool EdgeOrSideLength, typename Done>
567 bool
568 kamada_kawai_spring_layout(
569 const Graph& g,
570 PositionMap position,
571 WeightMap weight,
572 const Topology& topology,
573 detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length,
574 Done done,
575 typename property_traits<WeightMap>::value_type spring_constant)
576 {
577 return kamada_kawai_spring_layout(g, position, weight, topology, edge_or_side_length,
578 done, spring_constant,
579 get(vertex_index, g));
580 }
581
582 /**
583 * \overload
584 */
585 template<typename Topology, typename Graph, typename PositionMap, typename WeightMap,
586 typename T, bool EdgeOrSideLength, typename Done>
587 bool
588 kamada_kawai_spring_layout(
589 const Graph& g,
590 PositionMap position,
591 WeightMap weight,
592 const Topology& topology,
593 detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length,
594 Done done)
595 {
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));
599 }
600
601 /**
602 * \overload
603 */
604 template<typename Topology, typename Graph, typename PositionMap, typename WeightMap,
605 typename T, bool EdgeOrSideLength>
606 bool
607 kamada_kawai_spring_layout(
608 const Graph& g,
609 PositionMap position,
610 WeightMap weight,
611 const Topology& topology,
612 detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length)
613 {
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>(),
617 weight_type(1.0),
618 get(vertex_index, g));
619 }
620} // end namespace boost
621
622#endif // BOOST_GRAPH_KAMADA_KAWAI_SPRING_LAYOUT_HPP
623

source code of boost/boost/graph/kamada_kawai_spring_layout.hpp