#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
+namespace boost
+{
+namespace detail
+{
+ namespace graph
{
- 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;
+ /**
+ * 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;
}
- }
- 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[3][3], 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);
+ /**
+ * 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;
}
- return 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[3][3], 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 partial derivatives dE/dx_m and dE/dy_m
- deriv_type
- compute_partial_derivatives(vertex_descriptor m)
- {
+ // 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;
+ 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;
- }
+ 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;
- }
+ return result;
+ }
- // The actual Kamada-Kawai spring layout algorithm implementation
- bool run()
- {
+ // 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;
+ 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);
-
- // std::cerr << "edge_length = " << edge_length << std::endl;
-
- // Compute l_{ij} and k_{ij}
- const weight_type K = spring_constant;
- vertex_iterator ui, end;
- for (ui = vertices(g).first, end = vertices(g).second; 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);
+ deriv_type result;
- for (ui = vertices(g).first, end = vertices(g).second; ui != end; ++ui) {
- deriv_type deriv = compute_partial_derivatives(*ui);
- put(partial_derivatives, *ui, deriv);
+ // TBD: looks like an accumulate to me
+ BGL_FORALL_VERTICES_T(i, g, Graph)
+ {
+ deriv_type deriv = compute_partial_derivative(m, i);
+ result += deriv;
+ }
- weight_type delta = topology.norm(deriv);
+ return result;
+ }
- if (delta > delta_p) {
- p = *ui;
- delta_p = delta;
- }
- }
+ // The actual Kamada-Kawai spring layout algorithm implementation
+ bool run()
+ {
+#ifndef BOOST_NO_STDC_NAMESPACE
+ using std::sqrt;
+#endif // BOOST_NO_STDC_NAMESPACE
- 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, end = vertices(g).second; ui != end; ++ui) {
- vertex_descriptor i = *ui;
- p_partials[get(index, i)] = compute_partial_derivative(i, p);
- }
-
- do {
- // For debugging, compute the energy value E
- double E = 0.;
- for (ui = vertices(g).first, end = vertices(g).second; ui != end; ++ui) {
- vertex_iterator vi = ui;
- for (++vi; vi != end; ++vi) {
- double dist = topology.distance(position[*ui], position[*vi]);
- weight_type k_ij = spring_strength[get(index,*ui)][get(index,*vi)];
- weight_type l_ij = distance[get(index, *ui)][get(index, *vi)];
- E += .5 * k_ij * (dist - l_ij) * (dist - l_ij);
- }
- }
- // std::cerr << "E = " << E << std::endl;
-
- // 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, end = vertices(g).second; 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;
- // dE_d_d[i][j] += k_mi * l_mi * sqrt(hypot(diff[i], diff[j])) * inv_dist_cubed;
+ // 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);
+
+ // std::cerr << "edge_length = " << edge_length << std::endl;
+
+ // Compute l_{ij} and k_{ij}
+ const weight_type K = spring_constant;
+ vertex_iterator ui, end;
+ for (ui = vertices(g).first, end = vertices(g).second;
+ 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);
}
- }
}
- }
- }
-
- 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);
+ // Compute Delta_i and find max
+ vertex_descriptor p = *vertices(g).first;
+ weight_type delta_p(0);
- // Move p by delta
- position[p] = topology.adjust(position[p], delta);
+ for (ui = vertices(g).first, end = vertices(g).second;
+ ui != end; ++ui)
+ {
+ deriv_type deriv = compute_partial_derivatives(*ui);
+ put(partial_derivatives, *ui, deriv);
- // Recompute partial derivatives and delta_p
- deriv_type deriv = compute_partial_derivatives(p);
- put(partial_derivatives, p, deriv);
+ weight_type delta = topology.norm(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, end = vertices(g).second; 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);
+ if (delta > delta_p)
+ {
+ p = *ui;
+ delta_p = delta;
+ }
+ }
- deriv += old_p_partial - old_deriv_p;
+ 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, end = vertices(g).second;
+ ui != end; ++ui)
+ {
+ vertex_descriptor i = *ui;
+ p_partials[get(index, i)]
+ = compute_partial_derivative(i, p);
+ }
- put(partial_derivatives, *ui, deriv);
- weight_type delta = topology.norm(deriv);
+ do
+ {
+ // For debugging, compute the energy value E
+ double E = 0.;
+ for (ui = vertices(g).first, end = vertices(g).second;
+ ui != end; ++ui)
+ {
+ vertex_iterator vi = ui;
+ for (++vi; vi != end; ++vi)
+ {
+ double dist = topology.distance(
+ position[*ui], position[*vi]);
+ weight_type k_ij = spring_strength[get(
+ index, *ui)][get(index, *vi)];
+ weight_type l_ij = distance[get(index, *ui)]
+ [get(index, *vi)];
+ E += .5 * k_ij * (dist - l_ij) * (dist - l_ij);
+ }
+ }
+ // std::cerr << "E = " << E << std::endl;
+
+ // 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, end = vertices(g).second;
+ 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;
+ // dE_d_d[i][j] += k_mi * l_mi *
+ // sqrt(hypot(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, end = vertices(g).second;
+ 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;
+ }
+ }
+ }
- if (delta > delta_p) {
- p = *ui;
- delta_p = delta;
+ return true;
}
- }
- }
- 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
- {
+ 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)
+ : 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;
+ 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;
}
-
- 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);
+ 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;
}
-
- 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:
+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 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 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,
+};
+
+/** \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 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 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,
+ 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);
+{
+ 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,
+}
+
+/**
+ * \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,
+ typename property_traits< WeightMap >::value_type spring_constant,
VertexIndexMap index)
- {
- typedef typename property_traits<WeightMap>::value_type weight_type;
+{
+ 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;
+ 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);
+ 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(),
+ 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,
+ 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,
+ 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));
- }
+{
+ 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
projects / ceph.git / blobdiff