williamr@2: // Copyright 2004 The Trustees of Indiana University.
williamr@2: 
williamr@2: // Distributed under the Boost Software License, Version 1.0.
williamr@2: // (See accompanying file LICENSE_1_0.txt or copy at
williamr@2: // http://www.boost.org/LICENSE_1_0.txt)
williamr@2: 
williamr@2: //  Authors: Jeremiah Willcock
williamr@2: //           Douglas Gregor
williamr@2: //           Andrew Lumsdaine
williamr@2: #ifndef BOOST_GRAPH_GURSOY_ATUN_LAYOUT_HPP
williamr@2: #define BOOST_GRAPH_GURSOY_ATUN_LAYOUT_HPP
williamr@2: 
williamr@2: // Gursoy-Atun graph layout, based on:
williamr@2: // "Neighbourhood Preserving Load Balancing: A Self-Organizing Approach"
williamr@2: // in EuroPar 2000, p. 234 of LNCS 1900
williamr@2: // http://springerlink.metapress.com/link.asp?id=pcu07ew5rhexp9yt
williamr@2: 
williamr@2: #include <cmath>
williamr@2: #include <vector>
williamr@2: #include <exception>
williamr@2: #include <algorithm>
williamr@2: 
williamr@2: #include <boost/graph/visitors.hpp>
williamr@2: #include <boost/graph/properties.hpp>
williamr@2: #include <boost/random/uniform_01.hpp>
williamr@2: #include <boost/random/linear_congruential.hpp>
williamr@2: #include <boost/shared_ptr.hpp>
williamr@2: #include <boost/graph/breadth_first_search.hpp>
williamr@2: #include <boost/graph/dijkstra_shortest_paths.hpp>
williamr@2: #include <boost/graph/named_function_params.hpp>
williamr@2: 
williamr@2: namespace boost { 
williamr@2: 
williamr@2: namespace detail {
williamr@2: 
williamr@2: struct over_distance_limit : public std::exception {};
williamr@2: 
williamr@2: template <typename PositionMap, typename NodeDistanceMap,  typename Topology,
williamr@2:           typename Graph>
williamr@2: struct update_position_visitor {
williamr@2:   typedef typename Topology::point_type Point;
williamr@2:   PositionMap position_map;
williamr@2:   NodeDistanceMap node_distance;
williamr@2:   const Topology& space;
williamr@2:   Point input_vector;
williamr@2:   double distance_limit;
williamr@2:   double learning_constant;
williamr@2:   double falloff_ratio;
williamr@2: 
williamr@2:   typedef boost::on_examine_vertex event_filter;
williamr@2: 
williamr@2:   typedef typename graph_traits<Graph>::vertex_descriptor
williamr@2:     vertex_descriptor;
williamr@2: 
williamr@2:   update_position_visitor(PositionMap position_map,
williamr@2:                           NodeDistanceMap node_distance,
williamr@2:                           const Topology& space,
williamr@2:                           const Point& input_vector,
williamr@2:                           double distance_limit,
williamr@2:                           double learning_constant,
williamr@2:                           double falloff_ratio):
williamr@2:     position_map(position_map), node_distance(node_distance), 
williamr@2:     space(space),
williamr@2:     input_vector(input_vector), distance_limit(distance_limit),
williamr@2:     learning_constant(learning_constant), falloff_ratio(falloff_ratio) {}
williamr@2: 
williamr@2:   void operator()(vertex_descriptor v, const Graph&) const 
williamr@2:   {
williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE
williamr@2:     using std::pow;
williamr@2: #endif
williamr@2: 
williamr@2:     if (get(node_distance, v) > distance_limit)
williamr@2:       throw over_distance_limit();
williamr@2:     Point old_position = get(position_map, v);
williamr@2:     double distance = get(node_distance, v);
williamr@2:     double fraction = 
williamr@2:       learning_constant * pow(falloff_ratio, distance * distance);
williamr@2:     put(position_map, v,
williamr@2:         space.move_position_toward(old_position, fraction, input_vector));
williamr@2:   }
williamr@2: };
williamr@2: 
williamr@2: template<typename EdgeWeightMap>
williamr@2: struct gursoy_shortest
williamr@2: {
williamr@2:   template<typename Graph, typename NodeDistanceMap, typename UpdatePosition>
williamr@2:   static inline void 
williamr@2:   run(const Graph& g, typename graph_traits<Graph>::vertex_descriptor s,
williamr@2:       NodeDistanceMap node_distance,  UpdatePosition& update_position,
williamr@2:       EdgeWeightMap weight)
williamr@2:   {
williamr@2:     boost::dijkstra_shortest_paths(g, s, weight_map(weight).
williamr@2:       visitor(boost::make_dijkstra_visitor(std::make_pair(
williamr@2:        boost::record_distances(node_distance, boost::on_edge_relaxed()),
williamr@2:         update_position))));
williamr@2:   }
williamr@2: };
williamr@2: 
williamr@2: template<>
williamr@2: struct gursoy_shortest<dummy_property_map>
williamr@2: {
williamr@2:   template<typename Graph, typename NodeDistanceMap, typename UpdatePosition>
williamr@2:   static inline void 
williamr@2:   run(const Graph& g, typename graph_traits<Graph>::vertex_descriptor s,
williamr@2:       NodeDistanceMap node_distance,  UpdatePosition& update_position,
williamr@2:       dummy_property_map)
williamr@2:   {
williamr@2:     boost::breadth_first_search(g, s,
williamr@2:       visitor(boost::make_bfs_visitor(std::make_pair(
williamr@2:         boost::record_distances(node_distance, boost::on_tree_edge()),
williamr@2:         update_position))));
williamr@2:   }
williamr@2: };
williamr@2: 
williamr@2: } // namespace detail
williamr@2: 
williamr@2: template <typename VertexListAndIncidenceGraph,  typename Topology,
williamr@2:           typename PositionMap, typename Diameter, typename VertexIndexMap, 
williamr@2:           typename EdgeWeightMap>
williamr@2: void 
williamr@2: gursoy_atun_step
williamr@2:   (const VertexListAndIncidenceGraph& graph,  
williamr@2:    const Topology& space,
williamr@2:    PositionMap position,
williamr@2:    Diameter diameter,
williamr@2:    double learning_constant,
williamr@2:    VertexIndexMap vertex_index_map,
williamr@2:    EdgeWeightMap weight)
williamr@2: {
williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE
williamr@2:   using std::pow;
williamr@2:   using std::exp;
williamr@2: #endif
williamr@2: 
williamr@2:   typedef typename graph_traits<VertexListAndIncidenceGraph>::vertex_iterator
williamr@2:     vertex_iterator;
williamr@2:   typedef typename graph_traits<VertexListAndIncidenceGraph>::vertex_descriptor
williamr@2:     vertex_descriptor;
williamr@2:   typedef typename Topology::point_type point_type;
williamr@2:   vertex_iterator i, iend;
williamr@2:   std::vector<double> distance_from_input_vector(num_vertices(graph));
williamr@2:   typedef boost::iterator_property_map<std::vector<double>::iterator, 
williamr@2:                                        VertexIndexMap,
williamr@2:                                        double, double&>
williamr@2:     DistanceFromInputMap;
williamr@2:   DistanceFromInputMap distance_from_input(distance_from_input_vector.begin(),
williamr@2:                                            vertex_index_map);
williamr@2:   std::vector<double> node_distance_map_vector(num_vertices(graph));
williamr@2:   typedef boost::iterator_property_map<std::vector<double>::iterator, 
williamr@2:                                        VertexIndexMap,
williamr@2:                                        double, double&>
williamr@2:     NodeDistanceMap;
williamr@2:   NodeDistanceMap node_distance(node_distance_map_vector.begin(),
williamr@2:                                 vertex_index_map);
williamr@2:   point_type input_vector = space.random_point();
williamr@2:   vertex_descriptor min_distance_loc 
williamr@2:     = graph_traits<VertexListAndIncidenceGraph>::null_vertex();
williamr@2:   double min_distance = 0.0;
williamr@2:   bool min_distance_unset = true;
williamr@2:   for (boost::tie(i, iend) = vertices(graph); i != iend; ++i) {
williamr@2:     double this_distance = space.distance(get(position, *i), input_vector);
williamr@2:     put(distance_from_input, *i, this_distance);
williamr@2:     if (min_distance_unset || this_distance < min_distance) {
williamr@2:       min_distance = this_distance;
williamr@2:       min_distance_loc = *i;
williamr@2:     }
williamr@2:     min_distance_unset = false;
williamr@2:   }
williamr@2:   assert (!min_distance_unset); // Graph must have at least one vertex
williamr@2:   boost::detail::update_position_visitor<
williamr@2:       PositionMap, NodeDistanceMap, Topology,
williamr@2:       VertexListAndIncidenceGraph> 
williamr@2:     update_position(position, node_distance, space,
williamr@2:                     input_vector, diameter, learning_constant, 
williamr@2:                     exp(-1. / (2 * diameter * diameter)));
williamr@2:   std::fill(node_distance_map_vector.begin(), node_distance_map_vector.end(), 0);
williamr@2:   try {
williamr@2:     typedef detail::gursoy_shortest<EdgeWeightMap> shortest;
williamr@2:     shortest::run(graph, min_distance_loc, node_distance, update_position,
williamr@2:                   weight);    
williamr@2:   } catch (detail::over_distance_limit) { 
williamr@2:     /* Thrown to break out of BFS or Dijkstra early */ 
williamr@2:   }
williamr@2: }
williamr@2: 
williamr@2: template <typename VertexListAndIncidenceGraph,  typename Topology,
williamr@2:           typename PositionMap, typename VertexIndexMap, 
williamr@2:           typename EdgeWeightMap>
williamr@2: void gursoy_atun_refine(const VertexListAndIncidenceGraph& graph,  
williamr@2:                         const Topology& space,
williamr@2:                         PositionMap position,
williamr@2:                         int nsteps,
williamr@2:                         double diameter_initial,
williamr@2:                         double diameter_final,
williamr@2:                         double learning_constant_initial,
williamr@2:                         double learning_constant_final,
williamr@2:                         VertexIndexMap vertex_index_map,
williamr@2:                         EdgeWeightMap weight) 
williamr@2: {
williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE
williamr@2:   using std::pow;
williamr@2:   using std::exp;
williamr@2: #endif
williamr@2: 
williamr@2:   typedef typename graph_traits<VertexListAndIncidenceGraph>::vertex_iterator
williamr@2:     vertex_iterator;
williamr@2:   typedef typename graph_traits<VertexListAndIncidenceGraph>::vertex_descriptor
williamr@2:     vertex_descriptor;
williamr@2:   typedef typename Topology::point_type point_type;
williamr@2:   vertex_iterator i, iend;
williamr@2:   double diameter_ratio = (double)diameter_final / diameter_initial;
williamr@2:   double learning_constant_ratio = 
williamr@2:     learning_constant_final / learning_constant_initial;
williamr@2:   std::vector<double> distance_from_input_vector(num_vertices(graph));
williamr@2:   typedef boost::iterator_property_map<std::vector<double>::iterator, 
williamr@2:                                        VertexIndexMap,
williamr@2:                                        double, double&>
williamr@2:     DistanceFromInputMap;
williamr@2:   DistanceFromInputMap distance_from_input(distance_from_input_vector.begin(),
williamr@2:                                            vertex_index_map);
williamr@2:   std::vector<int> node_distance_map_vector(num_vertices(graph));
williamr@2:   typedef boost::iterator_property_map<std::vector<int>::iterator, 
williamr@2:                                        VertexIndexMap, double, double&>
williamr@2:     NodeDistanceMap;
williamr@2:   NodeDistanceMap node_distance(node_distance_map_vector.begin(),
williamr@2:                                 vertex_index_map);
williamr@2:   for (int round = 0; round < nsteps; ++round) {
williamr@2:     double part_done = (double)round / (nsteps - 1);
williamr@2:     int diameter = (int)(diameter_initial * pow(diameter_ratio, part_done));
williamr@2:     double learning_constant = 
williamr@2:       learning_constant_initial * pow(learning_constant_ratio, part_done);
williamr@2:     gursoy_atun_step(graph, space, position, diameter, learning_constant, 
williamr@2:                      vertex_index_map, weight);
williamr@2:   }
williamr@2: }
williamr@2: 
williamr@2: template <typename VertexListAndIncidenceGraph,  typename Topology,
williamr@2:           typename PositionMap, typename VertexIndexMap, 
williamr@2:           typename EdgeWeightMap>
williamr@2: void gursoy_atun_layout(const VertexListAndIncidenceGraph& graph,  
williamr@2:                         const Topology& space,
williamr@2:                         PositionMap position,
williamr@2:                         int nsteps,
williamr@2:                         double diameter_initial,
williamr@2:                         double diameter_final,
williamr@2:                         double learning_constant_initial,
williamr@2:                         double learning_constant_final,
williamr@2:                         VertexIndexMap vertex_index_map,
williamr@2:                         EdgeWeightMap weight)
williamr@2: {
williamr@2:   typedef typename graph_traits<VertexListAndIncidenceGraph>::vertex_iterator
williamr@2:     vertex_iterator;
williamr@2:   vertex_iterator i, iend;
williamr@2:   for (boost::tie(i, iend) = vertices(graph); i != iend; ++i) {
williamr@2:     put(position, *i, space.random_point());
williamr@2:   }
williamr@2:   gursoy_atun_refine(graph, space,
williamr@2:                      position, nsteps,
williamr@2:                      diameter_initial, diameter_final, 
williamr@2:                      learning_constant_initial, learning_constant_final,
williamr@2:                      vertex_index_map, weight);
williamr@2: }
williamr@2: 
williamr@2: template <typename VertexListAndIncidenceGraph,  typename Topology,
williamr@2:           typename PositionMap, typename VertexIndexMap>
williamr@2: void gursoy_atun_layout(const VertexListAndIncidenceGraph& graph,  
williamr@2:                         const Topology& space,
williamr@2:                         PositionMap position,
williamr@2:                         int nsteps,
williamr@2:                         double diameter_initial,
williamr@2:                         double diameter_final,
williamr@2:                         double learning_constant_initial,
williamr@2:                         double learning_constant_final,
williamr@2:                         VertexIndexMap vertex_index_map)
williamr@2: {
williamr@2:   gursoy_atun_layout(graph, space, position, nsteps, 
williamr@2:                      diameter_initial, diameter_final, 
williamr@2:                      learning_constant_initial, learning_constant_final, 
williamr@2:                      vertex_index_map, dummy_property_map());
williamr@2: }
williamr@2: 
williamr@2: template <typename VertexListAndIncidenceGraph, typename Topology,
williamr@2:           typename PositionMap>
williamr@2: void gursoy_atun_layout(const VertexListAndIncidenceGraph& graph,  
williamr@2:                         const Topology& space,
williamr@2:                         PositionMap position,
williamr@2:                         int nsteps,
williamr@2:                         double diameter_initial,
williamr@2:                         double diameter_final = 1.0,
williamr@2:                         double learning_constant_initial = 0.8,
williamr@2:                         double learning_constant_final = 0.2)
williamr@2: { 
williamr@2:   gursoy_atun_layout(graph, space, position, nsteps, diameter_initial,
williamr@2:                      diameter_final, learning_constant_initial,
williamr@2:                      learning_constant_final, get(vertex_index, graph)); 
williamr@2: }
williamr@2: 
williamr@2: template <typename VertexListAndIncidenceGraph, typename Topology,
williamr@2:           typename PositionMap>
williamr@2: void gursoy_atun_layout(const VertexListAndIncidenceGraph& graph,  
williamr@2:                         const Topology& space,
williamr@2:                         PositionMap position,
williamr@2:                         int nsteps)
williamr@2: {
williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE
williamr@2:   using std::sqrt;
williamr@2: #endif
williamr@2: 
williamr@2:   gursoy_atun_layout(graph, space, position, nsteps, 
williamr@2:                      sqrt((double)num_vertices(graph)));
williamr@2: }
williamr@2: 
williamr@2: template <typename VertexListAndIncidenceGraph, typename Topology,
williamr@2:           typename PositionMap>
williamr@2: void gursoy_atun_layout(const VertexListAndIncidenceGraph& graph,  
williamr@2:                         const Topology& space,
williamr@2:                         PositionMap position)
williamr@2: {
williamr@2:   gursoy_atun_layout(graph, space, position, num_vertices(graph));
williamr@2: }
williamr@2: 
williamr@2: template<typename VertexListAndIncidenceGraph, typename Topology,
williamr@2:          typename PositionMap, typename P, typename T, typename R>
williamr@2: void 
williamr@2: gursoy_atun_layout(const VertexListAndIncidenceGraph& graph,  
williamr@2:                    const Topology& space,
williamr@2:                    PositionMap position,
williamr@2:                    const bgl_named_params<P,T,R>& params)
williamr@2: {
williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE
williamr@2:   using std::sqrt;
williamr@2: #endif
williamr@2: 
williamr@2:   std::pair<double, double> diam(sqrt(double(num_vertices(graph))), 1.0);
williamr@2:   std::pair<double, double> learn(0.8, 0.2);
williamr@2:   gursoy_atun_layout(graph, space, position,
williamr@2:                      choose_param(get_param(params, iterations_t()),
williamr@2:                                   num_vertices(graph)),
williamr@2:                      choose_param(get_param(params, diameter_range_t()), 
williamr@2:                                   diam).first,
williamr@2:                      choose_param(get_param(params, diameter_range_t()), 
williamr@2:                                   diam).second,
williamr@2:                      choose_param(get_param(params, learning_constant_range_t()), 
williamr@2:                                   learn).first,
williamr@2:                      choose_param(get_param(params, learning_constant_range_t()), 
williamr@2:                                   learn).second,
williamr@2:                      choose_const_pmap(get_param(params, vertex_index), graph,
williamr@2:                                        vertex_index),
williamr@2:                      choose_param(get_param(params, edge_weight), 
williamr@2:                                   dummy_property_map()));
williamr@2: }
williamr@2: 
williamr@2: /***********************************************************
williamr@2:  * Topologies                                              *
williamr@2:  ***********************************************************/
williamr@2: template<std::size_t Dims>
williamr@2: class convex_topology 
williamr@2: {
williamr@2:   struct point 
williamr@2:   {
williamr@2:     point() { }
williamr@2:     double& operator[](std::size_t i) {return values[i];}
williamr@2:     const double& operator[](std::size_t i) const {return values[i];}
williamr@2: 
williamr@2:   private:
williamr@2:     double values[Dims];
williamr@2:   };
williamr@2: 
williamr@2:  public:
williamr@2:   typedef point point_type;
williamr@2: 
williamr@2:   double distance(point a, point b) const 
williamr@2:   {
williamr@2:     double dist = 0;
williamr@2:     for (std::size_t i = 0; i < Dims; ++i) {
williamr@2:       double diff = b[i] - a[i];
williamr@2:       dist += diff * diff;
williamr@2:     }
williamr@2:     // Exact properties of the distance are not important, as long as
williamr@2:     // < on what this returns matches real distances
williamr@2:     return dist;
williamr@2:   }
williamr@2: 
williamr@2:   point move_position_toward(point a, double fraction, point b) const 
williamr@2:   {
williamr@2:     point result;
williamr@2:     for (std::size_t i = 0; i < Dims; ++i)
williamr@2:       result[i] = a[i] + (b[i] - a[i]) * fraction;
williamr@2:     return result;
williamr@2:   }
williamr@2: };
williamr@2: 
williamr@2: template<std::size_t Dims,
williamr@2:          typename RandomNumberGenerator = minstd_rand>
williamr@2: class hypercube_topology : public convex_topology<Dims>
williamr@2: {
williamr@2:   typedef uniform_01<RandomNumberGenerator, double> rand_t;
williamr@2: 
williamr@2:  public:
williamr@2:   typedef typename convex_topology<Dims>::point_type point_type;
williamr@2: 
williamr@2:   explicit hypercube_topology(double scaling = 1.0) 
williamr@2:     : gen_ptr(new RandomNumberGenerator), rand(new rand_t(*gen_ptr)), 
williamr@2:       scaling(scaling) 
williamr@2:   { }
williamr@2: 
williamr@2:   hypercube_topology(RandomNumberGenerator& gen, double scaling = 1.0) 
williamr@2:     : gen_ptr(), rand(new rand_t(gen)), scaling(scaling) { }
williamr@2:                      
williamr@2:   point_type random_point() const 
williamr@2:   {
williamr@2:     point_type p;
williamr@2:     for (std::size_t i = 0; i < Dims; ++i)
williamr@2:       p[i] = (*rand)() * scaling;
williamr@2:     return p;
williamr@2:   }
williamr@2: 
williamr@2:  private:
williamr@2:   shared_ptr<RandomNumberGenerator> gen_ptr;
williamr@2:   shared_ptr<rand_t> rand;
williamr@2:   double scaling;
williamr@2: };
williamr@2: 
williamr@2: template<typename RandomNumberGenerator = minstd_rand>
williamr@2: class square_topology : public hypercube_topology<2, RandomNumberGenerator>
williamr@2: {
williamr@2:   typedef hypercube_topology<2, RandomNumberGenerator> inherited;
williamr@2: 
williamr@2:  public:
williamr@2:   explicit square_topology(double scaling = 1.0) : inherited(scaling) { }
williamr@2:   
williamr@2:   square_topology(RandomNumberGenerator& gen, double scaling = 1.0) 
williamr@2:     : inherited(gen, scaling) { }
williamr@2: };
williamr@2: 
williamr@2: template<typename RandomNumberGenerator = minstd_rand>
williamr@2: class cube_topology : public hypercube_topology<3, RandomNumberGenerator>
williamr@2: {
williamr@2:   typedef hypercube_topology<3, RandomNumberGenerator> inherited;
williamr@2: 
williamr@2:  public:
williamr@2:   explicit cube_topology(double scaling = 1.0) : inherited(scaling) { }
williamr@2:   
williamr@2:   cube_topology(RandomNumberGenerator& gen, double scaling = 1.0) 
williamr@2:     : inherited(gen, scaling) { }
williamr@2: };
williamr@2: 
williamr@2: template<std::size_t Dims,
williamr@2:          typename RandomNumberGenerator = minstd_rand>
williamr@2: class ball_topology : public convex_topology<Dims>
williamr@2: {
williamr@2:   typedef uniform_01<RandomNumberGenerator, double> rand_t;
williamr@2: 
williamr@2:  public:
williamr@2:   typedef typename convex_topology<Dims>::point_type point_type;
williamr@2: 
williamr@2:   explicit ball_topology(double radius = 1.0) 
williamr@2:     : gen_ptr(new RandomNumberGenerator), rand(new rand_t(*gen_ptr)), 
williamr@2:       radius(radius) 
williamr@2:   { }
williamr@2: 
williamr@2:   ball_topology(RandomNumberGenerator& gen, double radius = 1.0) 
williamr@2:     : gen_ptr(), rand(new rand_t(gen)), radius(radius) { }
williamr@2:                      
williamr@2:   point_type random_point() const 
williamr@2:   {
williamr@2:     point_type p;
williamr@2:     double dist_sum;
williamr@2:     do {
williamr@2:       dist_sum = 0.0;
williamr@2:       for (std::size_t i = 0; i < Dims; ++i) {
williamr@2:         double x = (*rand)() * 2*radius - radius;
williamr@2:         p[i] = x;
williamr@2:         dist_sum += x * x;
williamr@2:       }
williamr@2:     } while (dist_sum > radius*radius);
williamr@2:     return p;
williamr@2:   }
williamr@2: 
williamr@2:  private:
williamr@2:   shared_ptr<RandomNumberGenerator> gen_ptr;
williamr@2:   shared_ptr<rand_t> rand;
williamr@2:   double radius;
williamr@2: };
williamr@2: 
williamr@2: template<typename RandomNumberGenerator = minstd_rand>
williamr@2: class circle_topology : public ball_topology<2, RandomNumberGenerator>
williamr@2: {
williamr@2:   typedef ball_topology<2, RandomNumberGenerator> inherited;
williamr@2: 
williamr@2:  public:
williamr@2:   explicit circle_topology(double radius = 1.0) : inherited(radius) { }
williamr@2:   
williamr@2:   circle_topology(RandomNumberGenerator& gen, double radius = 1.0) 
williamr@2:     : inherited(gen, radius) { }
williamr@2: };
williamr@2: 
williamr@2: template<typename RandomNumberGenerator = minstd_rand>
williamr@2: class sphere_topology : public ball_topology<3, RandomNumberGenerator>
williamr@2: {
williamr@2:   typedef ball_topology<3, RandomNumberGenerator> inherited;
williamr@2: 
williamr@2:  public:
williamr@2:   explicit sphere_topology(double radius = 1.0) : inherited(radius) { }
williamr@2:   
williamr@2:   sphere_topology(RandomNumberGenerator& gen, double radius = 1.0) 
williamr@2:     : inherited(gen, radius) { }
williamr@2: };
williamr@2: 
williamr@2: template<typename RandomNumberGenerator = minstd_rand>
williamr@2: class heart_topology 
williamr@2: {
williamr@2:   // Heart is defined as the union of three shapes:
williamr@2:   // Square w/ corners (+-1000, -1000), (0, 0), (0, -2000)
williamr@2:   // Circle centered at (-500, -500) radius 500*sqrt(2)
williamr@2:   // Circle centered at (500, -500) radius 500*sqrt(2)
williamr@2:   // Bounding box (-1000, -2000) - (1000, 500*(sqrt(2) - 1))
williamr@2: 
williamr@2:   struct point 
williamr@2:   {
williamr@2:     point() { values[0] = 0.0; values[1] = 0.0; }
williamr@2:     point(double x, double y) { values[0] = x; values[1] = y; }
williamr@2: 
williamr@2:     double& operator[](std::size_t i)       { return values[i]; }
williamr@2:     double  operator[](std::size_t i) const { return values[i]; }
williamr@2: 
williamr@2:   private:
williamr@2:     double values[2];
williamr@2:   };
williamr@2: 
williamr@2:   bool in_heart(point p) const 
williamr@2:   {
williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE
williamr@2:     using std::abs;
williamr@2:     using std::pow;
williamr@2: #endif
williamr@2: 
williamr@2:     if (p[1] < abs(p[0]) - 2000) return false; // Bottom
williamr@2:     if (p[1] <= -1000) return true; // Diagonal of square
williamr@2:     if (pow(p[0] - -500, 2) + pow(p[1] - -500, 2) <= 500000)
williamr@2:       return true; // Left circle
williamr@2:     if (pow(p[0] - 500, 2) + pow(p[1] - -500, 2) <= 500000)
williamr@2:       return true; // Right circle
williamr@2:     return false;
williamr@2:   }
williamr@2: 
williamr@2:   bool segment_within_heart(point p1, point p2) const 
williamr@2:   {
williamr@2:     // Assumes that p1 and p2 are within the heart
williamr@2:     if ((p1[0] < 0) == (p2[0] < 0)) return true; // Same side of symmetry line
williamr@2:     if (p1[0] == p2[0]) return true; // Vertical
williamr@2:     double slope = (p2[1] - p1[1]) / (p2[0] - p1[0]);
williamr@2:     double intercept = p1[1] - p1[0] * slope;
williamr@2:     if (intercept > 0) return false; // Crosses between circles
williamr@2:     return true;
williamr@2:   }
williamr@2: 
williamr@2:   typedef uniform_01<RandomNumberGenerator, double> rand_t;
williamr@2: 
williamr@2:  public:
williamr@2:   typedef point point_type;
williamr@2: 
williamr@2:   heart_topology() 
williamr@2:     : gen_ptr(new RandomNumberGenerator), rand(new rand_t(*gen_ptr)) { }
williamr@2: 
williamr@2:   heart_topology(RandomNumberGenerator& gen) 
williamr@2:     : gen_ptr(), rand(new rand_t(gen)) { }
williamr@2: 
williamr@2:   point random_point() const 
williamr@2:   {
williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE
williamr@2:     using std::sqrt;
williamr@2: #endif
williamr@2: 
williamr@2:     point result;
williamr@2:     double sqrt2 = sqrt(2.);
williamr@2:     do {
williamr@2:       result[0] = (*rand)() * (1000 + 1000 * sqrt2) - (500 + 500 * sqrt2);
williamr@2:       result[1] = (*rand)() * (2000 + 500 * (sqrt2 - 1)) - 2000;
williamr@2:     } while (!in_heart(result));
williamr@2:     return result;
williamr@2:   }
williamr@2: 
williamr@2:   double distance(point a, point b) const 
williamr@2:   {
williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE
williamr@2:     using std::sqrt;
williamr@2: #endif
williamr@2:     if (segment_within_heart(a, b)) {
williamr@2:       // Straight line
williamr@2:       return sqrt((b[0] - a[0]) * (b[0] - a[0]) + (b[1] - a[1]) * (b[1] - a[1]));
williamr@2:     } else {
williamr@2:       // Straight line bending around (0, 0)
williamr@2:       return sqrt(a[0] * a[0] + a[1] * a[1]) + sqrt(b[0] * b[0] + b[1] * b[1]);
williamr@2:     }
williamr@2:   }
williamr@2: 
williamr@2:   point move_position_toward(point a, double fraction, point b) const 
williamr@2:   {
williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE
williamr@2:     using std::sqrt;
williamr@2: #endif
williamr@2: 
williamr@2:     if (segment_within_heart(a, b)) {
williamr@2:       // Straight line
williamr@2:       return point(a[0] + (b[0] - a[0]) * fraction,
williamr@2:                    a[1] + (b[1] - a[1]) * fraction);
williamr@2:     } else {
williamr@2:       double distance_to_point_a = sqrt(a[0] * a[0] + a[1] * a[1]);
williamr@2:       double distance_to_point_b = sqrt(b[0] * b[0] + b[1] * b[1]);
williamr@2:       double location_of_point = distance_to_point_a / 
williamr@2:                                    (distance_to_point_a + distance_to_point_b);
williamr@2:       if (fraction < location_of_point)
williamr@2:         return point(a[0] * (1 - fraction / location_of_point), 
williamr@2:                      a[1] * (1 - fraction / location_of_point));
williamr@2:       else
williamr@2:         return point(
williamr@2:           b[0] * ((fraction - location_of_point) / (1 - location_of_point)),
williamr@2:           b[1] * ((fraction - location_of_point) / (1 - location_of_point)));
williamr@2:     }
williamr@2:   }
williamr@2: 
williamr@2:  private:
williamr@2:   shared_ptr<RandomNumberGenerator> gen_ptr;
williamr@2:   shared_ptr<rand_t> rand;
williamr@2: };
williamr@2: 
williamr@2: } // namespace boost
williamr@2: 
williamr@2: #endif // BOOST_GRAPH_GURSOY_ATUN_LAYOUT_HPP