Chris@16: // Copyright 2004 The Trustees of Indiana University.
Chris@16: 
Chris@16: // Distributed under the Boost Software License, Version 1.0.
Chris@16: // (See accompanying file LICENSE_1_0.txt or copy at
Chris@16: // http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: //  Authors: Douglas Gregor
Chris@16: //           Andrew Lumsdaine
Chris@16: #ifndef BOOST_GRAPH_KAMADA_KAWAI_SPRING_LAYOUT_HPP
Chris@16: #define BOOST_GRAPH_KAMADA_KAWAI_SPRING_LAYOUT_HPP
Chris@16: 
Chris@16: #include <boost/graph/graph_traits.hpp>
Chris@16: #include <boost/graph/topology.hpp>
Chris@16: #include <boost/graph/iteration_macros.hpp>
Chris@16: #include <boost/graph/johnson_all_pairs_shortest.hpp>
Chris@16: #include <boost/type_traits/is_convertible.hpp>
Chris@16: #include <utility>
Chris@16: #include <iterator>
Chris@16: #include <vector>
Chris@16: #include <iostream>
Chris@16: #include <boost/limits.hpp>
Chris@16: #include <boost/config/no_tr1/cmath.hpp>
Chris@16: 
Chris@16: namespace boost {
Chris@16:   namespace detail { namespace graph {
Chris@16:     /**
Chris@16:      * Denotes an edge or display area side length used to scale a
Chris@16:      * Kamada-Kawai drawing.
Chris@16:      */
Chris@16:     template<bool Edge, typename T>
Chris@16:     struct edge_or_side
Chris@16:     {
Chris@16:       explicit edge_or_side(T value) : value(value) {}
Chris@16: 
Chris@16:       T value;
Chris@16:     };
Chris@16: 
Chris@16:     /**
Chris@16:      * Compute the edge length from an edge length. This is trivial.
Chris@16:      */
Chris@16:     template<typename Graph, typename DistanceMap, typename IndexMap, 
Chris@16:              typename T>
Chris@16:     T compute_edge_length(const Graph&, DistanceMap, IndexMap, 
Chris@16:                           edge_or_side<true, T> length)
Chris@16:     { return length.value; }
Chris@16: 
Chris@16:     /**
Chris@16:      * Compute the edge length based on the display area side
Chris@16:        length. We do this by dividing the side length by the largest
Chris@16:        shortest distance between any two vertices in the graph.
Chris@16:      */
Chris@16:     template<typename Graph, typename DistanceMap, typename IndexMap, 
Chris@16:              typename T>
Chris@16:     T
Chris@16:     compute_edge_length(const Graph& g, DistanceMap distance, IndexMap index,
Chris@16:                         edge_or_side<false, T> length)
Chris@16:     {
Chris@16:       T result(0);
Chris@16: 
Chris@16:       typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
Chris@16: 
Chris@16:       for (vertex_iterator ui = vertices(g).first, end = vertices(g).second;
Chris@16:            ui != end; ++ui) {
Chris@16:         vertex_iterator vi = ui;
Chris@16:         for (++vi; vi != end; ++vi) {
Chris@16:           T dij = distance[get(index, *ui)][get(index, *vi)];
Chris@16:           if (dij > result) result = dij;
Chris@16:         }
Chris@16:       }
Chris@16:       return length.value / result;
Chris@16:     }
Chris@16: 
Chris@16:     /**
Chris@16:      * Dense linear solver for fixed-size matrices.
Chris@16:      */
Chris@16:     template <std::size_t Size>
Chris@16:     struct linear_solver {
Chris@16:       // Indices in mat are (row, column)
Chris@16:       // template <typename Vec>
Chris@16:       // static Vec solve(double mat[Size][Size], Vec rhs);
Chris@16:     };
Chris@16: 
Chris@16:     template <>
Chris@16:     struct linear_solver<1> {
Chris@16:       template <typename Vec>
Chris@16:       static Vec solve(double mat[1][1], Vec rhs) {
Chris@16:         return rhs / mat[0][0];
Chris@16:       }
Chris@16:     };
Chris@16: 
Chris@16:     // These are from http://en.wikipedia.org/wiki/Cramer%27s_rule
Chris@16:     template <>
Chris@16:     struct linear_solver<2> {
Chris@16:       template <typename Vec>
Chris@16:       static Vec solve(double mat[2][2], Vec rhs) {
Chris@16:         double denom = mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1];
Chris@16:         double x_num = rhs[0]    * mat[1][1] - rhs[1]    * mat[0][1];
Chris@16:         double y_num = mat[0][0] * rhs[1]    - mat[1][0] * rhs[0]   ;
Chris@16:         Vec result;
Chris@16:         result[0] = x_num / denom;
Chris@16:         result[1] = y_num / denom;
Chris@16:         return result;
Chris@16:       }
Chris@16:     };
Chris@16: 
Chris@16:     template <>
Chris@16:     struct linear_solver<3> {
Chris@16:       template <typename Vec>
Chris@16:       static Vec solve(double mat[2][2], Vec rhs) {
Chris@16:         double denom = mat[0][0] * (mat[1][1] * mat[2][2] - mat[2][1] * mat[1][2])
Chris@16:                      - mat[1][0] * (mat[0][1] * mat[2][2] - mat[2][1] * mat[0][2])
Chris@16:                      + mat[2][0] * (mat[0][1] * mat[1][2] - mat[1][1] * mat[0][2]);
Chris@16:         double x_num = rhs[0]    * (mat[1][1] * mat[2][2] - mat[2][1] * mat[1][2])
Chris@16:                      - rhs[1]    * (mat[0][1] * mat[2][2] - mat[2][1] * mat[0][2])
Chris@16:                      + rhs[2]    * (mat[0][1] * mat[1][2] - mat[1][1] * mat[0][2]);
Chris@16:         double y_num = mat[0][0] * (rhs[1]    * mat[2][2] - rhs[2]    * mat[1][2])
Chris@16:                      - mat[1][0] * (rhs[0]    * mat[2][2] - rhs[2]    * mat[0][2])
Chris@16:                      + mat[2][0] * (rhs[0]    * mat[1][2] - rhs[1]    * mat[0][2]);
Chris@16:         double z_num = mat[0][0] * (mat[1][1] * rhs[2]    - mat[2][1] * rhs[1]   )
Chris@16:                      - mat[1][0] * (mat[0][1] * rhs[2]    - mat[2][1] * rhs[0]   )
Chris@16:                      + mat[2][0] * (mat[0][1] * rhs[1]    - mat[1][1] * rhs[0]   );
Chris@16:         Vec result;
Chris@16:         result[0] = x_num / denom;
Chris@16:         result[1] = y_num / denom;
Chris@16:         result[2] = z_num / denom;
Chris@16:         return result;
Chris@16:       }
Chris@16:     };
Chris@16: 
Chris@16:     /**
Chris@16:      * Implementation of the Kamada-Kawai spring layout algorithm.
Chris@16:      */
Chris@16:     template<typename Topology, typename Graph, typename PositionMap, typename WeightMap,
Chris@16:              typename EdgeOrSideLength, typename Done,
Chris@16:              typename VertexIndexMap, typename DistanceMatrix,
Chris@16:              typename SpringStrengthMatrix, typename PartialDerivativeMap>
Chris@16:     struct kamada_kawai_spring_layout_impl
Chris@16:     {
Chris@16:       typedef typename property_traits<WeightMap>::value_type weight_type;
Chris@16:       typedef typename Topology::point_type Point;
Chris@16:       typedef typename Topology::point_difference_type point_difference_type;
Chris@16:       typedef point_difference_type deriv_type;
Chris@16:       typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
Chris@16:       typedef typename graph_traits<Graph>::vertex_descriptor
Chris@16:         vertex_descriptor;
Chris@16: 
Chris@16:       kamada_kawai_spring_layout_impl(
Chris@16:         const Topology& topology,
Chris@16:         const Graph& g, 
Chris@16:         PositionMap position,
Chris@16:         WeightMap weight, 
Chris@16:         EdgeOrSideLength edge_or_side_length,
Chris@16:         Done done,
Chris@16:         weight_type spring_constant,
Chris@16:         VertexIndexMap index,
Chris@16:         DistanceMatrix distance,
Chris@16:         SpringStrengthMatrix spring_strength,
Chris@16:         PartialDerivativeMap partial_derivatives)
Chris@16:         : topology(topology), g(g), position(position), weight(weight), 
Chris@16:           edge_or_side_length(edge_or_side_length), done(done),
Chris@16:           spring_constant(spring_constant), index(index), distance(distance),
Chris@16:           spring_strength(spring_strength), 
Chris@16:           partial_derivatives(partial_derivatives) {}
Chris@16: 
Chris@16:       // Compute contribution of vertex i to the first partial
Chris@16:       // derivatives (dE/dx_m, dE/dy_m) (for vertex m)
Chris@16:       deriv_type
Chris@16:       compute_partial_derivative(vertex_descriptor m, vertex_descriptor i)
Chris@16:       {
Chris@16: #ifndef BOOST_NO_STDC_NAMESPACE
Chris@16:         using std::sqrt;
Chris@16: #endif // BOOST_NO_STDC_NAMESPACE
Chris@16: 
Chris@16:         deriv_type result;
Chris@16:         if (i != m) {
Chris@16:           point_difference_type diff = topology.difference(position[m], position[i]);
Chris@16:           weight_type dist = topology.norm(diff);
Chris@16:           result = spring_strength[get(index, m)][get(index, i)] 
Chris@16:             * (diff - distance[get(index, m)][get(index, i)]/dist*diff);
Chris@16:         }
Chris@16: 
Chris@16:         return result;
Chris@16:       }
Chris@16: 
Chris@16:       // Compute partial derivatives dE/dx_m and dE/dy_m
Chris@16:       deriv_type 
Chris@16:       compute_partial_derivatives(vertex_descriptor m)
Chris@16:       {
Chris@16: #ifndef BOOST_NO_STDC_NAMESPACE
Chris@16:         using std::sqrt;
Chris@16: #endif // BOOST_NO_STDC_NAMESPACE
Chris@16: 
Chris@16:         deriv_type result;
Chris@16: 
Chris@16:         // TBD: looks like an accumulate to me
Chris@16:         BGL_FORALL_VERTICES_T(i, g, Graph) {
Chris@16:           deriv_type deriv = compute_partial_derivative(m, i);
Chris@16:           result += deriv;
Chris@16:         }
Chris@16: 
Chris@16:         return result;
Chris@16:       }
Chris@16: 
Chris@16:       // The actual Kamada-Kawai spring layout algorithm implementation
Chris@16:       bool run()
Chris@16:       {
Chris@16: #ifndef BOOST_NO_STDC_NAMESPACE
Chris@16:         using std::sqrt;
Chris@16: #endif // BOOST_NO_STDC_NAMESPACE
Chris@16: 
Chris@16:         // Compute d_{ij} and place it in the distance matrix
Chris@16:         if (!johnson_all_pairs_shortest_paths(g, distance, index, weight, 
Chris@16:                                               weight_type(0)))
Chris@16:           return false;
Chris@16: 
Chris@16:         // Compute L based on side length (if needed), or retrieve L
Chris@16:         weight_type edge_length = 
Chris@16:           detail::graph::compute_edge_length(g, distance, index,
Chris@16:                                              edge_or_side_length);
Chris@16: 
Chris@16:         // std::cerr << "edge_length = " << edge_length << std::endl;
Chris@16:         
Chris@16:         // Compute l_{ij} and k_{ij}
Chris@16:         const weight_type K = spring_constant;
Chris@16:         vertex_iterator ui, end;
Chris@16:         for (ui = vertices(g).first, end = vertices(g).second; ui != end; ++ui) {
Chris@16:           vertex_iterator vi = ui;
Chris@16:           for (++vi; vi != end; ++vi) {
Chris@16:             weight_type dij = distance[get(index, *ui)][get(index, *vi)];
Chris@16:             if (dij == (std::numeric_limits<weight_type>::max)())
Chris@16:               return false;
Chris@16:             distance[get(index, *ui)][get(index, *vi)] = edge_length * dij;
Chris@16:             distance[get(index, *vi)][get(index, *ui)] = edge_length * dij;
Chris@16:             spring_strength[get(index, *ui)][get(index, *vi)] = K/(dij*dij);
Chris@16:             spring_strength[get(index, *vi)][get(index, *ui)] = K/(dij*dij);
Chris@16:           }
Chris@16:         }
Chris@16:         
Chris@16:         // Compute Delta_i and find max
Chris@16:         vertex_descriptor p = *vertices(g).first;
Chris@16:         weight_type delta_p(0);
Chris@16: 
Chris@16:         for (ui = vertices(g).first, end = vertices(g).second; ui != end; ++ui) {
Chris@16:           deriv_type deriv = compute_partial_derivatives(*ui);
Chris@16:           put(partial_derivatives, *ui, deriv);
Chris@16: 
Chris@16:           weight_type delta = topology.norm(deriv);
Chris@16: 
Chris@16:           if (delta > delta_p) {
Chris@16:             p = *ui;
Chris@16:             delta_p = delta;
Chris@16:           }
Chris@16:         }
Chris@16: 
Chris@16:         while (!done(delta_p, p, g, true)) {
Chris@16:           // The contribution p makes to the partial derivatives of
Chris@16:           // each vertex. Computing this (at O(n) cost) allows us to
Chris@16:           // update the delta_i values in O(n) time instead of O(n^2)
Chris@16:           // time.
Chris@16:           std::vector<deriv_type> p_partials(num_vertices(g));
Chris@16:           for (ui = vertices(g).first, end = vertices(g).second; ui != end; ++ui) {
Chris@16:             vertex_descriptor i = *ui;
Chris@16:             p_partials[get(index, i)] = compute_partial_derivative(i, p);
Chris@16:           }
Chris@16: 
Chris@16:           do {
Chris@16:             // For debugging, compute the energy value E
Chris@16:             double E = 0.;
Chris@16:             for (ui = vertices(g).first, end = vertices(g).second; ui != end; ++ui) {
Chris@16:               vertex_iterator vi = ui;
Chris@16:               for (++vi; vi != end; ++vi) {
Chris@16:                 double dist = topology.distance(position[*ui], position[*vi]);
Chris@16:                 weight_type k_ij = spring_strength[get(index,*ui)][get(index,*vi)];
Chris@16:                 weight_type l_ij = distance[get(index, *ui)][get(index, *vi)];
Chris@16:                 E += .5 * k_ij * (dist - l_ij) * (dist - l_ij);
Chris@16:               }
Chris@16:             }
Chris@16:             // std::cerr << "E = " << E << std::endl;
Chris@16: 
Chris@16:             // Compute the elements of the Jacobian
Chris@16:             // From
Chris@16:             // http://www.cs.panam.edu/~rfowler/papers/1994_kumar_fowler_A_Spring_UTPACSTR.pdf
Chris@16:             // with the bugs fixed in the off-diagonal case
Chris@16:             weight_type dE_d_d[Point::dimensions][Point::dimensions];
Chris@16:             for (std::size_t i = 0; i < Point::dimensions; ++i)
Chris@16:               for (std::size_t j = 0; j < Point::dimensions; ++j)
Chris@16:                 dE_d_d[i][j] = 0.;
Chris@16:             for (ui = vertices(g).first, end = vertices(g).second; ui != end; ++ui) {
Chris@16:               vertex_descriptor i = *ui;
Chris@16:               if (i != p) {
Chris@16:                 point_difference_type diff = topology.difference(position[p], position[i]);
Chris@16:                 weight_type dist = topology.norm(diff);
Chris@16:                 weight_type dist_squared = dist * dist;
Chris@16:                 weight_type inv_dist_cubed = 1. / (dist_squared * dist);
Chris@16:                 weight_type k_mi = spring_strength[get(index,p)][get(index,i)];
Chris@16:                 weight_type l_mi = distance[get(index, p)][get(index, i)];
Chris@16:                 for (std::size_t i = 0; i < Point::dimensions; ++i) {
Chris@16:                   for (std::size_t j = 0; j < Point::dimensions; ++j) {
Chris@16:                     if (i == j) {
Chris@16:                       dE_d_d[i][i] += k_mi * (1 + (l_mi * (diff[i] * diff[i] - dist_squared) * inv_dist_cubed));
Chris@16:                     } else {
Chris@16:                       dE_d_d[i][j] += k_mi * l_mi * diff[i] * diff[j] * inv_dist_cubed;
Chris@16:                       // dE_d_d[i][j] += k_mi * l_mi * sqrt(hypot(diff[i], diff[j])) * inv_dist_cubed;
Chris@16:                     }
Chris@16:                   }
Chris@16:                 }
Chris@16:               }
Chris@16:             }
Chris@16: 
Chris@16:             deriv_type dE_d = get(partial_derivatives, p);
Chris@16: 
Chris@16:             // Solve dE_d_d * delta = -dE_d to get delta
Chris@16:             point_difference_type delta = -linear_solver<Point::dimensions>::solve(dE_d_d, dE_d);
Chris@16: 
Chris@16:             // Move p by delta
Chris@16:             position[p] = topology.adjust(position[p], delta);
Chris@16: 
Chris@16:             // Recompute partial derivatives and delta_p
Chris@16:             deriv_type deriv = compute_partial_derivatives(p);
Chris@16:             put(partial_derivatives, p, deriv);
Chris@16: 
Chris@16:             delta_p = topology.norm(deriv);
Chris@16:           } while (!done(delta_p, p, g, false));
Chris@16: 
Chris@16:           // Select new p by updating each partial derivative and delta
Chris@16:           vertex_descriptor old_p = p;
Chris@16:           for (ui = vertices(g).first, end = vertices(g).second; ui != end; ++ui) {
Chris@16:             deriv_type old_deriv_p = p_partials[get(index, *ui)];
Chris@16:             deriv_type old_p_partial = 
Chris@16:               compute_partial_derivative(*ui, old_p);
Chris@16:             deriv_type deriv = get(partial_derivatives, *ui);
Chris@16: 
Chris@16:             deriv += old_p_partial - old_deriv_p;
Chris@16: 
Chris@16:             put(partial_derivatives, *ui, deriv);
Chris@16:             weight_type delta = topology.norm(deriv);
Chris@16: 
Chris@16:             if (delta > delta_p) {
Chris@16:               p = *ui;
Chris@16:               delta_p = delta;
Chris@16:             }
Chris@16:           }
Chris@16:         }
Chris@16: 
Chris@16:         return true;
Chris@16:       }
Chris@16: 
Chris@16:       const Topology& topology;
Chris@16:       const Graph& g; 
Chris@16:       PositionMap position;
Chris@16:       WeightMap weight; 
Chris@16:       EdgeOrSideLength edge_or_side_length;
Chris@16:       Done done;
Chris@16:       weight_type spring_constant;
Chris@16:       VertexIndexMap index;
Chris@16:       DistanceMatrix distance;
Chris@16:       SpringStrengthMatrix spring_strength;
Chris@16:       PartialDerivativeMap partial_derivatives;
Chris@16:     };
Chris@16:   } } // end namespace detail::graph
Chris@16: 
Chris@16:   /// States that the given quantity is an edge length.
Chris@16:   template<typename T> 
Chris@16:   inline detail::graph::edge_or_side<true, T>
Chris@16:   edge_length(T x) 
Chris@16:   { return detail::graph::edge_or_side<true, T>(x); }
Chris@16: 
Chris@16:   /// States that the given quantity is a display area side length.
Chris@16:   template<typename T> 
Chris@16:   inline detail::graph::edge_or_side<false, T>
Chris@16:   side_length(T x) 
Chris@16:   { return detail::graph::edge_or_side<false, T>(x); }
Chris@16: 
Chris@16:   /** 
Chris@16:    * \brief Determines when to terminate layout of a particular graph based
Chris@16:    * on a given relative tolerance. 
Chris@16:    */
Chris@16:   template<typename T = double>
Chris@16:   struct layout_tolerance
Chris@16:   {
Chris@16:     layout_tolerance(const T& tolerance = T(0.001))
Chris@16:       : tolerance(tolerance), last_energy((std::numeric_limits<T>::max)()),
Chris@16:         last_local_energy((std::numeric_limits<T>::max)()) { }
Chris@16: 
Chris@16:     template<typename Graph>
Chris@16:     bool 
Chris@16:     operator()(T delta_p, 
Chris@16:                typename boost::graph_traits<Graph>::vertex_descriptor p,
Chris@16:                const Graph& g,
Chris@16:                bool global)
Chris@16:     {
Chris@16:       if (global) {
Chris@16:         if (last_energy == (std::numeric_limits<T>::max)()) {
Chris@16:           last_energy = delta_p;
Chris@16:           return false;
Chris@16:         }
Chris@16:           
Chris@16:         T diff = last_energy - delta_p;
Chris@16:         if (diff < T(0)) diff = -diff;
Chris@16:         bool done = (delta_p == T(0) || diff / last_energy < tolerance);
Chris@16:         last_energy = delta_p;
Chris@16:         return done;
Chris@16:       } else {
Chris@16:         if (last_local_energy == (std::numeric_limits<T>::max)()) {
Chris@16:           last_local_energy = delta_p;
Chris@16:           return delta_p == T(0);
Chris@16:         }
Chris@16:           
Chris@16:         T diff = last_local_energy - delta_p;
Chris@16:         bool done = (delta_p == T(0) || (diff / last_local_energy) < tolerance);
Chris@16:         last_local_energy = delta_p;
Chris@16:         return done;
Chris@16:       }
Chris@16:     }
Chris@16: 
Chris@16:   private:
Chris@16:     T tolerance;
Chris@16:     T last_energy;
Chris@16:     T last_local_energy;
Chris@16:   };
Chris@16: 
Chris@16:   /** \brief Kamada-Kawai spring layout for undirected graphs.
Chris@16:    *
Chris@16:    * This algorithm performs graph layout (in two dimensions) for
Chris@16:    * connected, undirected graphs. It operates by relating the layout
Chris@16:    * of graphs to a dynamic spring system and minimizing the energy
Chris@16:    * within that system. The strength of a spring between two vertices
Chris@16:    * is inversely proportional to the square of the shortest distance
Chris@16:    * (in graph terms) between those two vertices. Essentially,
Chris@16:    * vertices that are closer in the graph-theoretic sense (i.e., by
Chris@16:    * following edges) will have stronger springs and will therefore be
Chris@16:    * placed closer together.
Chris@16:    *
Chris@16:    * Prior to invoking this algorithm, it is recommended that the
Chris@16:    * vertices be placed along the vertices of a regular n-sided
Chris@16:    * polygon.
Chris@16:    *
Chris@16:    * \param g (IN) must be a model of Vertex List Graph, Edge List
Chris@16:    * Graph, and Incidence Graph and must be undirected.
Chris@16:    *
Chris@16:    * \param position (OUT) must be a model of Lvalue Property Map,
Chris@16:    * where the value type is a class containing fields @c x and @c y
Chris@16:    * that will be set to the @c x and @c y coordinates of each vertex.
Chris@16:    *
Chris@16:    * \param weight (IN) must be a model of Readable Property Map,
Chris@16:    * which provides the weight of each edge in the graph @p g.
Chris@16:    *
Chris@16:    * \param topology (IN) must be a topology object (see topology.hpp),
Chris@16:    * which provides operations on points and differences between them.
Chris@16:    *
Chris@16:    * \param edge_or_side_length (IN) provides either the unit length
Chris@16:    * @c e of an edge in the layout or the length of a side @c s of the
Chris@16:    * display area, and must be either @c boost::edge_length(e) or @c
Chris@16:    * boost::side_length(s), respectively.
Chris@16:    *
Chris@16:    * \param done (IN) is a 4-argument function object that is passed
Chris@16:    * the current value of delta_p (i.e., the energy of vertex @p p),
Chris@16:    * the vertex @p p, the graph @p g, and a boolean flag indicating
Chris@16:    * whether @p delta_p is the maximum energy in the system (when @c
Chris@16:    * true) or the energy of the vertex being moved. Defaults to @c
Chris@16:    * layout_tolerance instantiated over the value type of the weight
Chris@16:    * map.
Chris@16:    *
Chris@16:    * \param spring_constant (IN) is the constant multiplied by each
Chris@16:    * spring's strength. Larger values create systems with more energy
Chris@16:    * that can take longer to stabilize; smaller values create systems
Chris@16:    * with less energy that stabilize quickly but do not necessarily
Chris@16:    * result in pleasing layouts. The default value is 1.
Chris@16:    *
Chris@16:    * \param index (IN) is a mapping from vertices to index values
Chris@16:    * between 0 and @c num_vertices(g). The default is @c
Chris@16:    * get(vertex_index,g).
Chris@16:    *
Chris@16:    * \param distance (UTIL/OUT) will be used to store the distance
Chris@16:    * from every vertex to every other vertex, which is computed in the
Chris@16:    * first stages of the algorithm. This value's type must be a model
Chris@16:    * of BasicMatrix with value type equal to the value type of the
Chris@16:    * weight map. The default is a vector of vectors.
Chris@16:    *
Chris@16:    * \param spring_strength (UTIL/OUT) will be used to store the
Chris@16:    * strength of the spring between every pair of vertices. This
Chris@16:    * value's type must be a model of BasicMatrix with value type equal
Chris@16:    * to the value type of the weight map. The default is a vector of
Chris@16:    * vectors.
Chris@16:    *
Chris@16:    * \param partial_derivatives (UTIL) will be used to store the
Chris@16:    * partial derivates of each vertex with respect to the @c x and @c
Chris@16:    * y coordinates. This must be a Read/Write Property Map whose value
Chris@16:    * type is a pair with both types equivalent to the value type of
Chris@16:    * the weight map. The default is an iterator property map.
Chris@16:    *
Chris@16:    * \returns @c true if layout was successful or @c false if a
Chris@16:    * negative weight cycle was detected.
Chris@16:    */
Chris@16:   template<typename Topology, typename Graph, typename PositionMap, typename WeightMap,
Chris@16:            typename T, bool EdgeOrSideLength, typename Done,
Chris@16:            typename VertexIndexMap, typename DistanceMatrix,
Chris@16:            typename SpringStrengthMatrix, typename PartialDerivativeMap>
Chris@16:   bool 
Chris@16:   kamada_kawai_spring_layout(
Chris@16:     const Graph& g, 
Chris@16:     PositionMap position,
Chris@16:     WeightMap weight, 
Chris@16:     const Topology& topology,
Chris@16:     detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length,
Chris@16:     Done done,
Chris@16:     typename property_traits<WeightMap>::value_type spring_constant,
Chris@16:     VertexIndexMap index,
Chris@16:     DistanceMatrix distance,
Chris@16:     SpringStrengthMatrix spring_strength,
Chris@16:     PartialDerivativeMap partial_derivatives)
Chris@16:   {
Chris@16:     BOOST_STATIC_ASSERT((is_convertible<
Chris@16:                            typename graph_traits<Graph>::directed_category*,
Chris@16:                            undirected_tag*
Chris@16:                          >::value));
Chris@16: 
Chris@16:     detail::graph::kamada_kawai_spring_layout_impl<
Chris@16:       Topology, Graph, PositionMap, WeightMap, 
Chris@16:       detail::graph::edge_or_side<EdgeOrSideLength, T>, Done, VertexIndexMap, 
Chris@16:       DistanceMatrix, SpringStrengthMatrix, PartialDerivativeMap>
Chris@16:       alg(topology, g, position, weight, edge_or_side_length, done, spring_constant,
Chris@16:           index, distance, spring_strength, partial_derivatives);
Chris@16:     return alg.run();
Chris@16:   }
Chris@16: 
Chris@16:   /**
Chris@16:    * \overload
Chris@16:    */
Chris@16:   template<typename Topology, typename Graph, typename PositionMap, typename WeightMap,
Chris@16:            typename T, bool EdgeOrSideLength, typename Done, 
Chris@16:            typename VertexIndexMap>
Chris@16:   bool 
Chris@16:   kamada_kawai_spring_layout(
Chris@16:     const Graph& g, 
Chris@16:     PositionMap position,
Chris@16:     WeightMap weight, 
Chris@16:     const Topology& topology,
Chris@16:     detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length,
Chris@16:     Done done,
Chris@16:     typename property_traits<WeightMap>::value_type spring_constant,
Chris@16:     VertexIndexMap index)
Chris@16:   {
Chris@16:     typedef typename property_traits<WeightMap>::value_type weight_type;
Chris@16: 
Chris@16:     typename graph_traits<Graph>::vertices_size_type n = num_vertices(g);
Chris@16:     typedef std::vector<weight_type> weight_vec;
Chris@16: 
Chris@16:     std::vector<weight_vec> distance(n, weight_vec(n));
Chris@16:     std::vector<weight_vec> spring_strength(n, weight_vec(n));
Chris@16:     std::vector<typename Topology::point_difference_type> partial_derivatives(n);
Chris@16: 
Chris@16:     return 
Chris@16:       kamada_kawai_spring_layout(
Chris@16:         g, position, weight, topology, edge_or_side_length, done, spring_constant, index,
Chris@16:         distance.begin(),
Chris@16:         spring_strength.begin(),
Chris@16:         make_iterator_property_map(partial_derivatives.begin(), index,
Chris@16:                                    typename Topology::point_difference_type()));
Chris@16:   }
Chris@16: 
Chris@16:   /**
Chris@16:    * \overload
Chris@16:    */
Chris@16:   template<typename Topology, typename Graph, typename PositionMap, typename WeightMap,
Chris@16:            typename T, bool EdgeOrSideLength, typename Done>
Chris@16:   bool 
Chris@16:   kamada_kawai_spring_layout(
Chris@16:     const Graph& g, 
Chris@16:     PositionMap position,
Chris@16:     WeightMap weight, 
Chris@16:     const Topology& topology,
Chris@16:     detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length,
Chris@16:     Done done,
Chris@16:     typename property_traits<WeightMap>::value_type spring_constant)
Chris@16:   {
Chris@16:     return kamada_kawai_spring_layout(g, position, weight, topology, edge_or_side_length,
Chris@16:                                       done, spring_constant, 
Chris@16:                                       get(vertex_index, g));
Chris@16:   }
Chris@16: 
Chris@16:   /**
Chris@16:    * \overload
Chris@16:    */
Chris@16:   template<typename Topology, typename Graph, typename PositionMap, typename WeightMap,
Chris@16:            typename T, bool EdgeOrSideLength, typename Done>
Chris@16:   bool 
Chris@16:   kamada_kawai_spring_layout(
Chris@16:     const Graph& g, 
Chris@16:     PositionMap position,
Chris@16:     WeightMap weight, 
Chris@16:     const Topology& topology,
Chris@16:     detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length,
Chris@16:     Done done)
Chris@16:   {
Chris@16:     typedef typename property_traits<WeightMap>::value_type weight_type;
Chris@16:     return kamada_kawai_spring_layout(g, position, weight, topology, edge_or_side_length,
Chris@16:                                       done, weight_type(1)); 
Chris@16:   }
Chris@16: 
Chris@16:   /**
Chris@16:    * \overload
Chris@16:    */
Chris@16:   template<typename Topology, typename Graph, typename PositionMap, typename WeightMap,
Chris@16:            typename T, bool EdgeOrSideLength>
Chris@16:   bool 
Chris@16:   kamada_kawai_spring_layout(
Chris@16:     const Graph& g, 
Chris@16:     PositionMap position,
Chris@16:     WeightMap weight, 
Chris@16:     const Topology& topology,
Chris@16:     detail::graph::edge_or_side<EdgeOrSideLength, T> edge_or_side_length)
Chris@16:   {
Chris@16:     typedef typename property_traits<WeightMap>::value_type weight_type;
Chris@16:     return kamada_kawai_spring_layout(g, position, weight, topology, edge_or_side_length,
Chris@16:                                       layout_tolerance<weight_type>(),
Chris@16:                                       weight_type(1.0), 
Chris@16:                                       get(vertex_index, g));
Chris@16:   }
Chris@16: } // end namespace boost
Chris@16: 
Chris@16: #endif // BOOST_GRAPH_KAMADA_KAWAI_SPRING_LAYOUT_HPP