source: NonGTP/Boost/boost/graph/kamada_kawai_spring_layout.hpp @ 857

Revision 857, 20.6 KB checked in by igarcia, 19 years ago (diff)
Line 
1// Copyright 2004 The Trustees of Indiana University.
2
3// Use, modification and distribution is subject to the Boost Software
4// License, Version 1.0. (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/johnson_all_pairs_shortest.hpp>
14#include <boost/type_traits/is_convertible.hpp>
15#include <utility>
16#include <iterator>
17#include <vector>
18#include <boost/limits.hpp>
19#include <cmath>
20
21namespace boost {
22  namespace detail { namespace graph {
23    /**
24     * Denotes an edge or display area side length used to scale a
25     * Kamada-Kawai drawing.
26     */
27    template<bool Edge, typename T>
28    struct edge_or_side
29    {
30      explicit edge_or_side(T value) : value(value) {}
31
32      T value;
33    };
34
35    /**
36     * Compute the edge length from an edge length. This is trivial.
37     */
38    template<typename Graph, typename DistanceMap, typename IndexMap,
39             typename T>
40    T compute_edge_length(const Graph&, DistanceMap, IndexMap,
41                          edge_or_side<true, T> length)
42    { return length.value; }
43
44    /**
45     * Compute the edge length based on the display area side
46       length. We do this by dividing the side length by the largest
47       shortest distance between any two vertices in the graph.
48     */
49    template<typename Graph, typename DistanceMap, typename IndexMap,
50             typename T>
51    T
52    compute_edge_length(const Graph& g, DistanceMap distance, IndexMap index,
53                        edge_or_side<false, T> length)
54    {
55      T result(0);
56
57      typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
58
59      for (vertex_iterator ui = vertices(g).first, end = vertices(g).second;
60           ui != end; ++ui) {
61        vertex_iterator vi = ui;
62        for (++vi; vi != end; ++vi) {
63          T dij = distance[get(index, *ui)][get(index, *vi)];
64          if (dij > result) result = dij;
65        }
66      }
67      return length.value / result;
68    }
69
70    /**
71     * Implementation of the Kamada-Kawai spring layout algorithm.
72     */
73    template<typename Graph, typename PositionMap, typename WeightMap,
74             typename EdgeOrSideLength, typename Done,
75             typename VertexIndexMap, typename DistanceMatrix,
76             typename SpringStrengthMatrix, typename PartialDerivativeMap>
77    struct kamada_kawai_spring_layout_impl
78    {
79      typedef typename property_traits<WeightMap>::value_type weight_type;
80      typedef std::pair<weight_type, weight_type> deriv_type;
81      typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
82      typedef typename graph_traits<Graph>::vertex_descriptor
83        vertex_descriptor;
84
85      kamada_kawai_spring_layout_impl(
86        const Graph& g,
87        PositionMap position,
88        WeightMap weight,
89        EdgeOrSideLength edge_or_side_length,
90        Done done,
91        weight_type spring_constant,
92        VertexIndexMap index,
93        DistanceMatrix distance,
94        SpringStrengthMatrix spring_strength,
95        PartialDerivativeMap partial_derivatives)
96        : g(g), position(position), weight(weight),
97          edge_or_side_length(edge_or_side_length), done(done),
98          spring_constant(spring_constant), index(index), distance(distance),
99          spring_strength(spring_strength),
100          partial_derivatives(partial_derivatives) {}
101
102      // Compute contribution of vertex i to the first partial
103      // derivatives (dE/dx_m, dE/dy_m) (for vertex m)
104      deriv_type
105      compute_partial_derivative(vertex_descriptor m, vertex_descriptor i)
106      {
107#ifndef BOOST_NO_STDC_NAMESPACE
108        using std::sqrt;
109#endif // BOOST_NO_STDC_NAMESPACE
110
111        deriv_type result(0, 0);
112        if (i != m) {
113          weight_type x_diff = position[m].x - position[i].x;
114          weight_type y_diff = position[m].y - position[i].y;
115          weight_type dist = sqrt(x_diff * x_diff + y_diff * y_diff);
116          result.first = spring_strength[get(index, m)][get(index, i)]
117            * (x_diff - distance[get(index, m)][get(index, i)]*x_diff/dist);
118          result.second = spring_strength[get(index, m)][get(index, i)]
119            * (y_diff - distance[get(index, m)][get(index, i)]*y_diff/dist);
120        }
121
122        return result;
123      }
124
125      // Compute partial derivatives dE/dx_m and dE/dy_m
126      deriv_type
127      compute_partial_derivatives(vertex_descriptor m)
128      {
129#ifndef BOOST_NO_STDC_NAMESPACE
130        using std::sqrt;
131#endif // BOOST_NO_STDC_NAMESPACE
132
133        deriv_type result(0, 0);
134
135        // TBD: looks like an accumulate to me
136        std::pair<vertex_iterator, vertex_iterator> verts = vertices(g);
137        for (/* no init */; verts.first != verts.second; ++verts.first) {
138          vertex_descriptor i = *verts.first;
139          deriv_type deriv = compute_partial_derivative(m, i);
140          result.first += deriv.first;
141          result.second += deriv.second;
142        }
143
144        return result;
145      }
146
147      // The actual Kamada-Kawai spring layout algorithm implementation
148      bool run()
149      {
150#ifndef BOOST_NO_STDC_NAMESPACE
151        using std::sqrt;
152#endif // BOOST_NO_STDC_NAMESPACE
153
154        // Compute d_{ij} and place it in the distance matrix
155        if (!johnson_all_pairs_shortest_paths(g, distance, index, weight,
156                                              weight_type(0)))
157          return false;
158
159        // Compute L based on side length (if needed), or retrieve L
160        weight_type edge_length =
161          detail::graph::compute_edge_length(g, distance, index,
162                                             edge_or_side_length);
163       
164        // Compute l_{ij} and k_{ij}
165        const weight_type K = spring_constant;
166        vertex_iterator ui, end = vertices(g).second;
167        for (ui = vertices(g).first; ui != end; ++ui) {
168          vertex_iterator vi = ui;
169          for (++vi; vi != end; ++vi) {
170            weight_type dij = distance[get(index, *ui)][get(index, *vi)];
171            if (dij == (std::numeric_limits<weight_type>::max)())
172              return false;
173            distance[get(index, *ui)][get(index, *vi)] = edge_length * dij;
174            distance[get(index, *vi)][get(index, *ui)] = edge_length * dij;
175            spring_strength[get(index, *ui)][get(index, *vi)] = K/(dij*dij);
176            spring_strength[get(index, *vi)][get(index, *ui)] = K/(dij*dij);
177          }
178        }
179       
180        // Compute Delta_i and find max
181        vertex_descriptor p = *vertices(g).first;
182        weight_type delta_p(0);
183
184        for (ui = vertices(g).first; ui != end; ++ui) {
185          deriv_type deriv = compute_partial_derivatives(*ui);
186          put(partial_derivatives, *ui, deriv);
187
188          weight_type delta =
189            sqrt(deriv.first*deriv.first + deriv.second*deriv.second);
190
191          if (delta > delta_p) {
192            p = *ui;
193            delta_p = delta;
194          }
195        }
196
197        while (!done(delta_p, p, g, true)) {
198          // The contribution p makes to the partial derivatives of
199          // each vertex. Computing this (at O(n) cost) allows us to
200          // update the delta_i values in O(n) time instead of O(n^2)
201          // time.
202          std::vector<deriv_type> p_partials(num_vertices(g));
203          for (ui = vertices(g).first; ui != end; ++ui) {
204            vertex_descriptor i = *ui;
205            p_partials[get(index, i)] = compute_partial_derivative(i, p);
206          }
207
208          do {
209            // Compute the 4 elements of the Jacobian
210            weight_type dE_dx_dx = 0, dE_dx_dy = 0, dE_dy_dx = 0, dE_dy_dy = 0;
211            for (ui = vertices(g).first; ui != end; ++ui) {
212              vertex_descriptor i = *ui;
213              if (i != p) {
214                weight_type x_diff = position[p].x - position[i].x;
215                weight_type y_diff = position[p].y - position[i].y;
216                weight_type dist = sqrt(x_diff * x_diff + y_diff * y_diff);
217                weight_type dist_cubed = dist * dist * dist;
218                weight_type k_mi = spring_strength[get(index,p)][get(index,i)];
219                weight_type l_mi = distance[get(index, p)][get(index, i)];
220                dE_dx_dx += k_mi * (1 - (l_mi * y_diff * y_diff)/dist_cubed);
221                dE_dx_dy += k_mi * l_mi * x_diff * y_diff / dist_cubed;
222                dE_dy_dx += k_mi * l_mi * x_diff * y_diff / dist_cubed;
223                dE_dy_dy += k_mi * (1 - (l_mi * x_diff * x_diff)/dist_cubed);
224              }
225            }
226
227            // Solve for delta_x and delta_y
228            weight_type dE_dx = get(partial_derivatives, p).first;
229            weight_type dE_dy = get(partial_derivatives, p).second;
230
231            weight_type delta_x =
232              (dE_dx_dy * dE_dy - dE_dy_dy * dE_dx)
233              / (dE_dx_dx * dE_dy_dy - dE_dx_dy * dE_dy_dx);
234
235            weight_type delta_y =
236              (dE_dx_dx * dE_dy - dE_dy_dx * dE_dx)
237              / (dE_dy_dx * dE_dx_dy - dE_dx_dx * dE_dy_dy);
238
239
240            // Move p by (delta_x, delta_y)
241            position[p].x += delta_x;
242            position[p].y += delta_y;
243
244            // Recompute partial derivatives and delta_p
245            deriv_type deriv = compute_partial_derivatives(p);
246            put(partial_derivatives, p, deriv);
247
248            delta_p =
249              sqrt(deriv.first*deriv.first + deriv.second*deriv.second);
250          } while (!done(delta_p, p, g, false));
251
252          // Select new p by updating each partial derivative and delta
253          vertex_descriptor old_p = p;
254          for (ui = vertices(g).first; ui != end; ++ui) {
255            deriv_type old_deriv_p = p_partials[get(index, *ui)];
256            deriv_type old_p_partial =
257              compute_partial_derivative(*ui, old_p);
258            deriv_type deriv = get(partial_derivatives, *ui);
259
260            deriv.first += old_p_partial.first - old_deriv_p.first;
261            deriv.second += old_p_partial.second - old_deriv_p.second;
262
263            put(partial_derivatives, *ui, deriv);
264            weight_type delta =
265              sqrt(deriv.first*deriv.first + deriv.second*deriv.second);
266
267            if (delta > delta_p) {
268              p = *ui;
269              delta_p = delta;
270            }
271          }
272        }
273
274        return true;
275      }
276
277      const Graph& g;
278      PositionMap position;
279      WeightMap weight;
280      EdgeOrSideLength edge_or_side_length;
281      Done done;
282      weight_type spring_constant;
283      VertexIndexMap index;
284      DistanceMatrix distance;
285      SpringStrengthMatrix spring_strength;
286      PartialDerivativeMap partial_derivatives;
287    };
288  } } // end namespace detail::graph
289
290  /// States that the given quantity is an edge length.
291  template<typename T>
292  inline detail::graph::edge_or_side<true, T>
293  edge_length(T x)
294  { return detail::graph::edge_or_side<true, T>(x); }
295
296  /// States that the given quantity is a display area side length.
297  template<typename T>
298  inline detail::graph::edge_or_side<false, T>
299  side_length(T x)
300  { return detail::graph::edge_or_side<false, T>(x); }
301
302  /**
303   * \brief Determines when to terminate layout of a particular graph based
304   * on a given relative tolerance.
305   */
306  template<typename T = double>
307  struct layout_tolerance
308  {
309    layout_tolerance(const T& tolerance = T(0.001))
310      : tolerance(tolerance), last_energy((std::numeric_limits<T>::max)()),
311        last_local_energy((std::numeric_limits<T>::max)()) { }
312
313    template<typename Graph>
314    bool
315    operator()(T delta_p,
316               typename boost::graph_traits<Graph>::vertex_descriptor p,
317               const Graph& g,
318               bool global)
319    {
320      if (global) {
321        if (last_energy == (std::numeric_limits<T>::max)()) {
322          last_energy = delta_p;
323          return false;
324        }
325         
326        T diff = last_energy - delta_p;
327        if (diff < T(0)) diff = -diff;
328        bool done = (delta_p == T(0) || diff / last_energy < tolerance);
329        last_energy = delta_p;
330        return done;
331      } else {
332        if (last_local_energy == (std::numeric_limits<T>::max)()) {
333          last_local_energy = delta_p;
334          return delta_p == T(0);
335        }
336         
337        T diff = last_local_energy - delta_p;
338        bool done = (delta_p == T(0) || (diff / last_local_energy) < tolerance);
339        last_local_energy = delta_p;
340        return done;
341      }
342    }
343
344  private:
345    T tolerance;
346    T last_energy;
347    T last_local_energy;
348  };
349
350  /** \brief Kamada-Kawai spring layout for undirected graphs.
351   *
352   * This algorithm performs graph layout (in two dimensions) for
353   * connected, undirected graphs. It operates by relating the layout
354   * of graphs to a dynamic spring system and minimizing the energy
355   * within that system. The strength of a spring between two vertices
356   * is inversely proportional to the square of the shortest distance
357   * (in graph terms) between those two vertices. Essentially,
358   * vertices that are closer in the graph-theoretic sense (i.e., by
359   * following edges) will have stronger springs and will therefore be
360   * placed closer together.
361   *
362   * Prior to invoking this algorithm, it is recommended that the
363   * vertices be placed along the vertices of a regular n-sided
364   * polygon.
365   *
366   * \param g (IN) must be a model of Vertex List Graph, Edge List
367   * Graph, and Incidence Graph and must be undirected.
368   *
369   * \param position (OUT) must be a model of Lvalue Property Map,
370   * where the value type is a class containing fields @c x and @c y
371   * that will be set to the @c x and @c y coordinates of each vertex.
372   *
373   * \param weight (IN) must be a model of Readable Property Map,
374   * which provides the weight of each edge in the graph @p g.
375   *
376   * \param edge_or_side_length (IN) provides either the unit length
377   * @c e of an edge in the layout or the length of a side @c s of the
378   * display area, and must be either @c boost::edge_length(e) or @c
379   * boost::side_length(s), respectively.
380   *
381   * \param done (IN) is a 4-argument function object that is passed
382   * the current value of delta_p (i.e., the energy of vertex @p p),
383   * the vertex @p p, the graph @p g, and a boolean flag indicating
384   * whether @p delta_p is the maximum energy in the system (when @c
385   * true) or the energy of the vertex being moved. Defaults to @c
386   * layout_tolerance instantiated over the value type of the weight
387   * map.
388   *
389   * \param spring_constant (IN) is the constant multiplied by each
390   * spring's strength. Larger values create systems with more energy
391   * that can take longer to stabilize; smaller values create systems
392   * with less energy that stabilize quickly but do not necessarily
393   * result in pleasing layouts. The default value is 1.
394   *
395   * \param index (IN) is a mapping from vertices to index values
396   * between 0 and @c num_vertices(g). The default is @c
397   * get(vertex_index,g).
398   *
399   * \param distance (UTIL/OUT) will be used to store the distance
400   * from every vertex to every other vertex, which is computed in the
401   * first stages of the algorithm. This value's type must be a model
402   * of BasicMatrix with value type equal to the value type of the
403   * weight map. The default is a a vector of vectors.
404   *
405   * \param spring_strength (UTIL/OUT) will be used to store the
406   * strength of the spring between every pair of vertices. This
407   * value's type must be a model of BasicMatrix with value type equal
408   * to the value type of the weight map. The default is a a vector of
409   * vectors.
410   *
411   * \param partial_derivatives (UTIL) will be used to store the
412   * partial derivates of each vertex with respect to the @c x and @c
413   * y coordinates. This must be a Read/Write Property Map whose value
414   * type is a pair with both types equivalent to the value type of
415   * the weight map. The default is an iterator property map.
416   *
417   * \returns @c true if layout was successful or @c false if a
418   * negative weight cycle was detected.
419   */
420  template<typename Graph, typename PositionMap, typename WeightMap,
421           typename T, bool EdgeOrSideLength, typename Done,
422           typename VertexIndexMap, typename DistanceMatrix,
423           typename SpringStrengthMatrix, typename PartialDerivativeMap>
424  bool
425  kamada_kawai_spring_layout(
426    const Graph& g,
427    PositionMap position,
428    WeightMap weight,
429    detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length,
430    Done done,
431    typename property_traits<WeightMap>::value_type spring_constant,
432    VertexIndexMap index,
433    DistanceMatrix distance,
434    SpringStrengthMatrix spring_strength,
435    PartialDerivativeMap partial_derivatives)
436  {
437    BOOST_STATIC_ASSERT((is_convertible<
438                           typename graph_traits<Graph>::directed_category*,
439                           undirected_tag*
440                         >::value));
441
442    detail::graph::kamada_kawai_spring_layout_impl<
443      Graph, PositionMap, WeightMap,
444      detail::graph::edge_or_side<EdgeOrSideLength, T>, Done, VertexIndexMap,
445      DistanceMatrix, SpringStrengthMatrix, PartialDerivativeMap>
446      alg(g, position, weight, edge_or_side_length, done, spring_constant,
447          index, distance, spring_strength, partial_derivatives);
448    return alg.run();
449  }
450
451  /**
452   * \overload
453   */
454  template<typename Graph, typename PositionMap, typename WeightMap,
455           typename T, bool EdgeOrSideLength, typename Done,
456           typename VertexIndexMap>
457  bool
458  kamada_kawai_spring_layout(
459    const Graph& g,
460    PositionMap position,
461    WeightMap weight,
462    detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length,
463    Done done,
464    typename property_traits<WeightMap>::value_type spring_constant,
465    VertexIndexMap index)
466  {
467    typedef typename property_traits<WeightMap>::value_type weight_type;
468
469    typename graph_traits<Graph>::vertices_size_type n = num_vertices(g);
470    typedef std::vector<weight_type> weight_vec;
471
472    std::vector<weight_vec> distance(n, weight_vec(n));
473    std::vector<weight_vec> spring_strength(n, weight_vec(n));
474    std::vector<std::pair<weight_type, weight_type> > partial_derivatives(n);
475
476    return
477      kamada_kawai_spring_layout(
478        g, position, weight, edge_or_side_length, done, spring_constant, index,
479        distance.begin(),
480        spring_strength.begin(),
481        make_iterator_property_map(partial_derivatives.begin(), index,
482                                   std::pair<weight_type, weight_type>()));
483  }
484
485  /**
486   * \overload
487   */
488  template<typename Graph, typename PositionMap, typename WeightMap,
489           typename T, bool EdgeOrSideLength, typename Done>
490  bool
491  kamada_kawai_spring_layout(
492    const Graph& g,
493    PositionMap position,
494    WeightMap weight,
495    detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length,
496    Done done,
497    typename property_traits<WeightMap>::value_type spring_constant)
498  {
499    return kamada_kawai_spring_layout(g, position, weight, edge_or_side_length,
500                                      done, spring_constant,
501                                      get(vertex_index, g));
502  }
503
504  /**
505   * \overload
506   */
507  template<typename Graph, typename PositionMap, typename WeightMap,
508           typename T, bool EdgeOrSideLength, typename Done>
509  bool
510  kamada_kawai_spring_layout(
511    const Graph& g,
512    PositionMap position,
513    WeightMap weight,
514    detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length,
515    Done done)
516  {
517    typedef typename property_traits<WeightMap>::value_type weight_type;
518    return kamada_kawai_spring_layout(g, position, weight, edge_or_side_length,
519                                      done, weight_type(1));
520  }
521
522  /**
523   * \overload
524   */
525  template<typename Graph, typename PositionMap, typename WeightMap,
526           typename T, bool EdgeOrSideLength>
527  bool
528  kamada_kawai_spring_layout(
529    const Graph& g,
530    PositionMap position,
531    WeightMap weight,
532    detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length)
533  {
534    typedef typename property_traits<WeightMap>::value_type weight_type;
535    return kamada_kawai_spring_layout(g, position, weight, edge_or_side_length,
536                                      layout_tolerance<weight_type>(),
537                                      weight_type(1.0),
538                                      get(vertex_index, g));
539  }
540} // end namespace boost
541
542#endif // BOOST_GRAPH_KAMADA_KAWAI_SPRING_LAYOUT_HPP
Note: See TracBrowser for help on using the repository browser.