ossrv_pub/boost_apis/boost/graph/gursoy_atun_layout.hpp
changeset 0 e4d67989cc36
equal deleted inserted replaced
-1:000000000000 0:e4d67989cc36
       
     1 // Copyright 2004 The Trustees of Indiana University.
       
     2 
       
     3 // Distributed under the Boost Software License, Version 1.0.
       
     4 // (See accompanying file LICENSE_1_0.txt or copy at
       
     5 // http://www.boost.org/LICENSE_1_0.txt)
       
     6 
       
     7 //  Authors: Jeremiah Willcock
       
     8 //           Douglas Gregor
       
     9 //           Andrew Lumsdaine
       
    10 #ifndef BOOST_GRAPH_GURSOY_ATUN_LAYOUT_HPP
       
    11 #define BOOST_GRAPH_GURSOY_ATUN_LAYOUT_HPP
       
    12 
       
    13 // Gursoy-Atun graph layout, based on:
       
    14 // "Neighbourhood Preserving Load Balancing: A Self-Organizing Approach"
       
    15 // in EuroPar 2000, p. 234 of LNCS 1900
       
    16 // http://springerlink.metapress.com/link.asp?id=pcu07ew5rhexp9yt
       
    17 
       
    18 #include <cmath>
       
    19 #include <vector>
       
    20 #include <exception>
       
    21 #include <algorithm>
       
    22 
       
    23 #include <boost/graph/visitors.hpp>
       
    24 #include <boost/graph/properties.hpp>
       
    25 #include <boost/random/uniform_01.hpp>
       
    26 #include <boost/random/linear_congruential.hpp>
       
    27 #include <boost/shared_ptr.hpp>
       
    28 #include <boost/graph/breadth_first_search.hpp>
       
    29 #include <boost/graph/dijkstra_shortest_paths.hpp>
       
    30 #include <boost/graph/named_function_params.hpp>
       
    31 
       
    32 namespace boost { 
       
    33 
       
    34 namespace detail {
       
    35 
       
    36 struct over_distance_limit : public std::exception {};
       
    37 
       
    38 template <typename PositionMap, typename NodeDistanceMap,  typename Topology,
       
    39           typename Graph>
       
    40 struct update_position_visitor {
       
    41   typedef typename Topology::point_type Point;
       
    42   PositionMap position_map;
       
    43   NodeDistanceMap node_distance;
       
    44   const Topology& space;
       
    45   Point input_vector;
       
    46   double distance_limit;
       
    47   double learning_constant;
       
    48   double falloff_ratio;
       
    49 
       
    50   typedef boost::on_examine_vertex event_filter;
       
    51 
       
    52   typedef typename graph_traits<Graph>::vertex_descriptor
       
    53     vertex_descriptor;
       
    54 
       
    55   update_position_visitor(PositionMap position_map,
       
    56                           NodeDistanceMap node_distance,
       
    57                           const Topology& space,
       
    58                           const Point& input_vector,
       
    59                           double distance_limit,
       
    60                           double learning_constant,
       
    61                           double falloff_ratio):
       
    62     position_map(position_map), node_distance(node_distance), 
       
    63     space(space),
       
    64     input_vector(input_vector), distance_limit(distance_limit),
       
    65     learning_constant(learning_constant), falloff_ratio(falloff_ratio) {}
       
    66 
       
    67   void operator()(vertex_descriptor v, const Graph&) const 
       
    68   {
       
    69 #ifndef BOOST_NO_STDC_NAMESPACE
       
    70     using std::pow;
       
    71 #endif
       
    72 
       
    73     if (get(node_distance, v) > distance_limit)
       
    74       throw over_distance_limit();
       
    75     Point old_position = get(position_map, v);
       
    76     double distance = get(node_distance, v);
       
    77     double fraction = 
       
    78       learning_constant * pow(falloff_ratio, distance * distance);
       
    79     put(position_map, v,
       
    80         space.move_position_toward(old_position, fraction, input_vector));
       
    81   }
       
    82 };
       
    83 
       
    84 template<typename EdgeWeightMap>
       
    85 struct gursoy_shortest
       
    86 {
       
    87   template<typename Graph, typename NodeDistanceMap, typename UpdatePosition>
       
    88   static inline void 
       
    89   run(const Graph& g, typename graph_traits<Graph>::vertex_descriptor s,
       
    90       NodeDistanceMap node_distance,  UpdatePosition& update_position,
       
    91       EdgeWeightMap weight)
       
    92   {
       
    93     boost::dijkstra_shortest_paths(g, s, weight_map(weight).
       
    94       visitor(boost::make_dijkstra_visitor(std::make_pair(
       
    95        boost::record_distances(node_distance, boost::on_edge_relaxed()),
       
    96         update_position))));
       
    97   }
       
    98 };
       
    99 
       
   100 template<>
       
   101 struct gursoy_shortest<dummy_property_map>
       
   102 {
       
   103   template<typename Graph, typename NodeDistanceMap, typename UpdatePosition>
       
   104   static inline void 
       
   105   run(const Graph& g, typename graph_traits<Graph>::vertex_descriptor s,
       
   106       NodeDistanceMap node_distance,  UpdatePosition& update_position,
       
   107       dummy_property_map)
       
   108   {
       
   109     boost::breadth_first_search(g, s,
       
   110       visitor(boost::make_bfs_visitor(std::make_pair(
       
   111         boost::record_distances(node_distance, boost::on_tree_edge()),
       
   112         update_position))));
       
   113   }
       
   114 };
       
   115 
       
   116 } // namespace detail
       
   117 
       
   118 template <typename VertexListAndIncidenceGraph,  typename Topology,
       
   119           typename PositionMap, typename Diameter, typename VertexIndexMap, 
       
   120           typename EdgeWeightMap>
       
   121 void 
       
   122 gursoy_atun_step
       
   123   (const VertexListAndIncidenceGraph& graph,  
       
   124    const Topology& space,
       
   125    PositionMap position,
       
   126    Diameter diameter,
       
   127    double learning_constant,
       
   128    VertexIndexMap vertex_index_map,
       
   129    EdgeWeightMap weight)
       
   130 {
       
   131 #ifndef BOOST_NO_STDC_NAMESPACE
       
   132   using std::pow;
       
   133   using std::exp;
       
   134 #endif
       
   135 
       
   136   typedef typename graph_traits<VertexListAndIncidenceGraph>::vertex_iterator
       
   137     vertex_iterator;
       
   138   typedef typename graph_traits<VertexListAndIncidenceGraph>::vertex_descriptor
       
   139     vertex_descriptor;
       
   140   typedef typename Topology::point_type point_type;
       
   141   vertex_iterator i, iend;
       
   142   std::vector<double> distance_from_input_vector(num_vertices(graph));
       
   143   typedef boost::iterator_property_map<std::vector<double>::iterator, 
       
   144                                        VertexIndexMap,
       
   145                                        double, double&>
       
   146     DistanceFromInputMap;
       
   147   DistanceFromInputMap distance_from_input(distance_from_input_vector.begin(),
       
   148                                            vertex_index_map);
       
   149   std::vector<double> node_distance_map_vector(num_vertices(graph));
       
   150   typedef boost::iterator_property_map<std::vector<double>::iterator, 
       
   151                                        VertexIndexMap,
       
   152                                        double, double&>
       
   153     NodeDistanceMap;
       
   154   NodeDistanceMap node_distance(node_distance_map_vector.begin(),
       
   155                                 vertex_index_map);
       
   156   point_type input_vector = space.random_point();
       
   157   vertex_descriptor min_distance_loc 
       
   158     = graph_traits<VertexListAndIncidenceGraph>::null_vertex();
       
   159   double min_distance = 0.0;
       
   160   bool min_distance_unset = true;
       
   161   for (boost::tie(i, iend) = vertices(graph); i != iend; ++i) {
       
   162     double this_distance = space.distance(get(position, *i), input_vector);
       
   163     put(distance_from_input, *i, this_distance);
       
   164     if (min_distance_unset || this_distance < min_distance) {
       
   165       min_distance = this_distance;
       
   166       min_distance_loc = *i;
       
   167     }
       
   168     min_distance_unset = false;
       
   169   }
       
   170   assert (!min_distance_unset); // Graph must have at least one vertex
       
   171   boost::detail::update_position_visitor<
       
   172       PositionMap, NodeDistanceMap, Topology,
       
   173       VertexListAndIncidenceGraph> 
       
   174     update_position(position, node_distance, space,
       
   175                     input_vector, diameter, learning_constant, 
       
   176                     exp(-1. / (2 * diameter * diameter)));
       
   177   std::fill(node_distance_map_vector.begin(), node_distance_map_vector.end(), 0);
       
   178   try {
       
   179     typedef detail::gursoy_shortest<EdgeWeightMap> shortest;
       
   180     shortest::run(graph, min_distance_loc, node_distance, update_position,
       
   181                   weight);    
       
   182   } catch (detail::over_distance_limit) { 
       
   183     /* Thrown to break out of BFS or Dijkstra early */ 
       
   184   }
       
   185 }
       
   186 
       
   187 template <typename VertexListAndIncidenceGraph,  typename Topology,
       
   188           typename PositionMap, typename VertexIndexMap, 
       
   189           typename EdgeWeightMap>
       
   190 void gursoy_atun_refine(const VertexListAndIncidenceGraph& graph,  
       
   191                         const Topology& space,
       
   192                         PositionMap position,
       
   193                         int nsteps,
       
   194                         double diameter_initial,
       
   195                         double diameter_final,
       
   196                         double learning_constant_initial,
       
   197                         double learning_constant_final,
       
   198                         VertexIndexMap vertex_index_map,
       
   199                         EdgeWeightMap weight) 
       
   200 {
       
   201 #ifndef BOOST_NO_STDC_NAMESPACE
       
   202   using std::pow;
       
   203   using std::exp;
       
   204 #endif
       
   205 
       
   206   typedef typename graph_traits<VertexListAndIncidenceGraph>::vertex_iterator
       
   207     vertex_iterator;
       
   208   typedef typename graph_traits<VertexListAndIncidenceGraph>::vertex_descriptor
       
   209     vertex_descriptor;
       
   210   typedef typename Topology::point_type point_type;
       
   211   vertex_iterator i, iend;
       
   212   double diameter_ratio = (double)diameter_final / diameter_initial;
       
   213   double learning_constant_ratio = 
       
   214     learning_constant_final / learning_constant_initial;
       
   215   std::vector<double> distance_from_input_vector(num_vertices(graph));
       
   216   typedef boost::iterator_property_map<std::vector<double>::iterator, 
       
   217                                        VertexIndexMap,
       
   218                                        double, double&>
       
   219     DistanceFromInputMap;
       
   220   DistanceFromInputMap distance_from_input(distance_from_input_vector.begin(),
       
   221                                            vertex_index_map);
       
   222   std::vector<int> node_distance_map_vector(num_vertices(graph));
       
   223   typedef boost::iterator_property_map<std::vector<int>::iterator, 
       
   224                                        VertexIndexMap, double, double&>
       
   225     NodeDistanceMap;
       
   226   NodeDistanceMap node_distance(node_distance_map_vector.begin(),
       
   227                                 vertex_index_map);
       
   228   for (int round = 0; round < nsteps; ++round) {
       
   229     double part_done = (double)round / (nsteps - 1);
       
   230     int diameter = (int)(diameter_initial * pow(diameter_ratio, part_done));
       
   231     double learning_constant = 
       
   232       learning_constant_initial * pow(learning_constant_ratio, part_done);
       
   233     gursoy_atun_step(graph, space, position, diameter, learning_constant, 
       
   234                      vertex_index_map, weight);
       
   235   }
       
   236 }
       
   237 
       
   238 template <typename VertexListAndIncidenceGraph,  typename Topology,
       
   239           typename PositionMap, typename VertexIndexMap, 
       
   240           typename EdgeWeightMap>
       
   241 void gursoy_atun_layout(const VertexListAndIncidenceGraph& graph,  
       
   242                         const Topology& space,
       
   243                         PositionMap position,
       
   244                         int nsteps,
       
   245                         double diameter_initial,
       
   246                         double diameter_final,
       
   247                         double learning_constant_initial,
       
   248                         double learning_constant_final,
       
   249                         VertexIndexMap vertex_index_map,
       
   250                         EdgeWeightMap weight)
       
   251 {
       
   252   typedef typename graph_traits<VertexListAndIncidenceGraph>::vertex_iterator
       
   253     vertex_iterator;
       
   254   vertex_iterator i, iend;
       
   255   for (boost::tie(i, iend) = vertices(graph); i != iend; ++i) {
       
   256     put(position, *i, space.random_point());
       
   257   }
       
   258   gursoy_atun_refine(graph, space,
       
   259                      position, nsteps,
       
   260                      diameter_initial, diameter_final, 
       
   261                      learning_constant_initial, learning_constant_final,
       
   262                      vertex_index_map, weight);
       
   263 }
       
   264 
       
   265 template <typename VertexListAndIncidenceGraph,  typename Topology,
       
   266           typename PositionMap, typename VertexIndexMap>
       
   267 void gursoy_atun_layout(const VertexListAndIncidenceGraph& graph,  
       
   268                         const Topology& space,
       
   269                         PositionMap position,
       
   270                         int nsteps,
       
   271                         double diameter_initial,
       
   272                         double diameter_final,
       
   273                         double learning_constant_initial,
       
   274                         double learning_constant_final,
       
   275                         VertexIndexMap vertex_index_map)
       
   276 {
       
   277   gursoy_atun_layout(graph, space, position, nsteps, 
       
   278                      diameter_initial, diameter_final, 
       
   279                      learning_constant_initial, learning_constant_final, 
       
   280                      vertex_index_map, dummy_property_map());
       
   281 }
       
   282 
       
   283 template <typename VertexListAndIncidenceGraph, typename Topology,
       
   284           typename PositionMap>
       
   285 void gursoy_atun_layout(const VertexListAndIncidenceGraph& graph,  
       
   286                         const Topology& space,
       
   287                         PositionMap position,
       
   288                         int nsteps,
       
   289                         double diameter_initial,
       
   290                         double diameter_final = 1.0,
       
   291                         double learning_constant_initial = 0.8,
       
   292                         double learning_constant_final = 0.2)
       
   293 { 
       
   294   gursoy_atun_layout(graph, space, position, nsteps, diameter_initial,
       
   295                      diameter_final, learning_constant_initial,
       
   296                      learning_constant_final, get(vertex_index, graph)); 
       
   297 }
       
   298 
       
   299 template <typename VertexListAndIncidenceGraph, typename Topology,
       
   300           typename PositionMap>
       
   301 void gursoy_atun_layout(const VertexListAndIncidenceGraph& graph,  
       
   302                         const Topology& space,
       
   303                         PositionMap position,
       
   304                         int nsteps)
       
   305 {
       
   306 #ifndef BOOST_NO_STDC_NAMESPACE
       
   307   using std::sqrt;
       
   308 #endif
       
   309 
       
   310   gursoy_atun_layout(graph, space, position, nsteps, 
       
   311                      sqrt((double)num_vertices(graph)));
       
   312 }
       
   313 
       
   314 template <typename VertexListAndIncidenceGraph, typename Topology,
       
   315           typename PositionMap>
       
   316 void gursoy_atun_layout(const VertexListAndIncidenceGraph& graph,  
       
   317                         const Topology& space,
       
   318                         PositionMap position)
       
   319 {
       
   320   gursoy_atun_layout(graph, space, position, num_vertices(graph));
       
   321 }
       
   322 
       
   323 template<typename VertexListAndIncidenceGraph, typename Topology,
       
   324          typename PositionMap, typename P, typename T, typename R>
       
   325 void 
       
   326 gursoy_atun_layout(const VertexListAndIncidenceGraph& graph,  
       
   327                    const Topology& space,
       
   328                    PositionMap position,
       
   329                    const bgl_named_params<P,T,R>& params)
       
   330 {
       
   331 #ifndef BOOST_NO_STDC_NAMESPACE
       
   332   using std::sqrt;
       
   333 #endif
       
   334 
       
   335   std::pair<double, double> diam(sqrt(double(num_vertices(graph))), 1.0);
       
   336   std::pair<double, double> learn(0.8, 0.2);
       
   337   gursoy_atun_layout(graph, space, position,
       
   338                      choose_param(get_param(params, iterations_t()),
       
   339                                   num_vertices(graph)),
       
   340                      choose_param(get_param(params, diameter_range_t()), 
       
   341                                   diam).first,
       
   342                      choose_param(get_param(params, diameter_range_t()), 
       
   343                                   diam).second,
       
   344                      choose_param(get_param(params, learning_constant_range_t()), 
       
   345                                   learn).first,
       
   346                      choose_param(get_param(params, learning_constant_range_t()), 
       
   347                                   learn).second,
       
   348                      choose_const_pmap(get_param(params, vertex_index), graph,
       
   349                                        vertex_index),
       
   350                      choose_param(get_param(params, edge_weight), 
       
   351                                   dummy_property_map()));
       
   352 }
       
   353 
       
   354 /***********************************************************
       
   355  * Topologies                                              *
       
   356  ***********************************************************/
       
   357 template<std::size_t Dims>
       
   358 class convex_topology 
       
   359 {
       
   360   struct point 
       
   361   {
       
   362     point() { }
       
   363     double& operator[](std::size_t i) {return values[i];}
       
   364     const double& operator[](std::size_t i) const {return values[i];}
       
   365 
       
   366   private:
       
   367     double values[Dims];
       
   368   };
       
   369 
       
   370  public:
       
   371   typedef point point_type;
       
   372 
       
   373   double distance(point a, point b) const 
       
   374   {
       
   375     double dist = 0;
       
   376     for (std::size_t i = 0; i < Dims; ++i) {
       
   377       double diff = b[i] - a[i];
       
   378       dist += diff * diff;
       
   379     }
       
   380     // Exact properties of the distance are not important, as long as
       
   381     // < on what this returns matches real distances
       
   382     return dist;
       
   383   }
       
   384 
       
   385   point move_position_toward(point a, double fraction, point b) const 
       
   386   {
       
   387     point result;
       
   388     for (std::size_t i = 0; i < Dims; ++i)
       
   389       result[i] = a[i] + (b[i] - a[i]) * fraction;
       
   390     return result;
       
   391   }
       
   392 };
       
   393 
       
   394 template<std::size_t Dims,
       
   395          typename RandomNumberGenerator = minstd_rand>
       
   396 class hypercube_topology : public convex_topology<Dims>
       
   397 {
       
   398   typedef uniform_01<RandomNumberGenerator, double> rand_t;
       
   399 
       
   400  public:
       
   401   typedef typename convex_topology<Dims>::point_type point_type;
       
   402 
       
   403   explicit hypercube_topology(double scaling = 1.0) 
       
   404     : gen_ptr(new RandomNumberGenerator), rand(new rand_t(*gen_ptr)), 
       
   405       scaling(scaling) 
       
   406   { }
       
   407 
       
   408   hypercube_topology(RandomNumberGenerator& gen, double scaling = 1.0) 
       
   409     : gen_ptr(), rand(new rand_t(gen)), scaling(scaling) { }
       
   410                      
       
   411   point_type random_point() const 
       
   412   {
       
   413     point_type p;
       
   414     for (std::size_t i = 0; i < Dims; ++i)
       
   415       p[i] = (*rand)() * scaling;
       
   416     return p;
       
   417   }
       
   418 
       
   419  private:
       
   420   shared_ptr<RandomNumberGenerator> gen_ptr;
       
   421   shared_ptr<rand_t> rand;
       
   422   double scaling;
       
   423 };
       
   424 
       
   425 template<typename RandomNumberGenerator = minstd_rand>
       
   426 class square_topology : public hypercube_topology<2, RandomNumberGenerator>
       
   427 {
       
   428   typedef hypercube_topology<2, RandomNumberGenerator> inherited;
       
   429 
       
   430  public:
       
   431   explicit square_topology(double scaling = 1.0) : inherited(scaling) { }
       
   432   
       
   433   square_topology(RandomNumberGenerator& gen, double scaling = 1.0) 
       
   434     : inherited(gen, scaling) { }
       
   435 };
       
   436 
       
   437 template<typename RandomNumberGenerator = minstd_rand>
       
   438 class cube_topology : public hypercube_topology<3, RandomNumberGenerator>
       
   439 {
       
   440   typedef hypercube_topology<3, RandomNumberGenerator> inherited;
       
   441 
       
   442  public:
       
   443   explicit cube_topology(double scaling = 1.0) : inherited(scaling) { }
       
   444   
       
   445   cube_topology(RandomNumberGenerator& gen, double scaling = 1.0) 
       
   446     : inherited(gen, scaling) { }
       
   447 };
       
   448 
       
   449 template<std::size_t Dims,
       
   450          typename RandomNumberGenerator = minstd_rand>
       
   451 class ball_topology : public convex_topology<Dims>
       
   452 {
       
   453   typedef uniform_01<RandomNumberGenerator, double> rand_t;
       
   454 
       
   455  public:
       
   456   typedef typename convex_topology<Dims>::point_type point_type;
       
   457 
       
   458   explicit ball_topology(double radius = 1.0) 
       
   459     : gen_ptr(new RandomNumberGenerator), rand(new rand_t(*gen_ptr)), 
       
   460       radius(radius) 
       
   461   { }
       
   462 
       
   463   ball_topology(RandomNumberGenerator& gen, double radius = 1.0) 
       
   464     : gen_ptr(), rand(new rand_t(gen)), radius(radius) { }
       
   465                      
       
   466   point_type random_point() const 
       
   467   {
       
   468     point_type p;
       
   469     double dist_sum;
       
   470     do {
       
   471       dist_sum = 0.0;
       
   472       for (std::size_t i = 0; i < Dims; ++i) {
       
   473         double x = (*rand)() * 2*radius - radius;
       
   474         p[i] = x;
       
   475         dist_sum += x * x;
       
   476       }
       
   477     } while (dist_sum > radius*radius);
       
   478     return p;
       
   479   }
       
   480 
       
   481  private:
       
   482   shared_ptr<RandomNumberGenerator> gen_ptr;
       
   483   shared_ptr<rand_t> rand;
       
   484   double radius;
       
   485 };
       
   486 
       
   487 template<typename RandomNumberGenerator = minstd_rand>
       
   488 class circle_topology : public ball_topology<2, RandomNumberGenerator>
       
   489 {
       
   490   typedef ball_topology<2, RandomNumberGenerator> inherited;
       
   491 
       
   492  public:
       
   493   explicit circle_topology(double radius = 1.0) : inherited(radius) { }
       
   494   
       
   495   circle_topology(RandomNumberGenerator& gen, double radius = 1.0) 
       
   496     : inherited(gen, radius) { }
       
   497 };
       
   498 
       
   499 template<typename RandomNumberGenerator = minstd_rand>
       
   500 class sphere_topology : public ball_topology<3, RandomNumberGenerator>
       
   501 {
       
   502   typedef ball_topology<3, RandomNumberGenerator> inherited;
       
   503 
       
   504  public:
       
   505   explicit sphere_topology(double radius = 1.0) : inherited(radius) { }
       
   506   
       
   507   sphere_topology(RandomNumberGenerator& gen, double radius = 1.0) 
       
   508     : inherited(gen, radius) { }
       
   509 };
       
   510 
       
   511 template<typename RandomNumberGenerator = minstd_rand>
       
   512 class heart_topology 
       
   513 {
       
   514   // Heart is defined as the union of three shapes:
       
   515   // Square w/ corners (+-1000, -1000), (0, 0), (0, -2000)
       
   516   // Circle centered at (-500, -500) radius 500*sqrt(2)
       
   517   // Circle centered at (500, -500) radius 500*sqrt(2)
       
   518   // Bounding box (-1000, -2000) - (1000, 500*(sqrt(2) - 1))
       
   519 
       
   520   struct point 
       
   521   {
       
   522     point() { values[0] = 0.0; values[1] = 0.0; }
       
   523     point(double x, double y) { values[0] = x; values[1] = y; }
       
   524 
       
   525     double& operator[](std::size_t i)       { return values[i]; }
       
   526     double  operator[](std::size_t i) const { return values[i]; }
       
   527 
       
   528   private:
       
   529     double values[2];
       
   530   };
       
   531 
       
   532   bool in_heart(point p) const 
       
   533   {
       
   534 #ifndef BOOST_NO_STDC_NAMESPACE
       
   535     using std::abs;
       
   536     using std::pow;
       
   537 #endif
       
   538 
       
   539     if (p[1] < abs(p[0]) - 2000) return false; // Bottom
       
   540     if (p[1] <= -1000) return true; // Diagonal of square
       
   541     if (pow(p[0] - -500, 2) + pow(p[1] - -500, 2) <= 500000)
       
   542       return true; // Left circle
       
   543     if (pow(p[0] - 500, 2) + pow(p[1] - -500, 2) <= 500000)
       
   544       return true; // Right circle
       
   545     return false;
       
   546   }
       
   547 
       
   548   bool segment_within_heart(point p1, point p2) const 
       
   549   {
       
   550     // Assumes that p1 and p2 are within the heart
       
   551     if ((p1[0] < 0) == (p2[0] < 0)) return true; // Same side of symmetry line
       
   552     if (p1[0] == p2[0]) return true; // Vertical
       
   553     double slope = (p2[1] - p1[1]) / (p2[0] - p1[0]);
       
   554     double intercept = p1[1] - p1[0] * slope;
       
   555     if (intercept > 0) return false; // Crosses between circles
       
   556     return true;
       
   557   }
       
   558 
       
   559   typedef uniform_01<RandomNumberGenerator, double> rand_t;
       
   560 
       
   561  public:
       
   562   typedef point point_type;
       
   563 
       
   564   heart_topology() 
       
   565     : gen_ptr(new RandomNumberGenerator), rand(new rand_t(*gen_ptr)) { }
       
   566 
       
   567   heart_topology(RandomNumberGenerator& gen) 
       
   568     : gen_ptr(), rand(new rand_t(gen)) { }
       
   569 
       
   570   point random_point() const 
       
   571   {
       
   572 #ifndef BOOST_NO_STDC_NAMESPACE
       
   573     using std::sqrt;
       
   574 #endif
       
   575 
       
   576     point result;
       
   577     double sqrt2 = sqrt(2.);
       
   578     do {
       
   579       result[0] = (*rand)() * (1000 + 1000 * sqrt2) - (500 + 500 * sqrt2);
       
   580       result[1] = (*rand)() * (2000 + 500 * (sqrt2 - 1)) - 2000;
       
   581     } while (!in_heart(result));
       
   582     return result;
       
   583   }
       
   584 
       
   585   double distance(point a, point b) const 
       
   586   {
       
   587 #ifndef BOOST_NO_STDC_NAMESPACE
       
   588     using std::sqrt;
       
   589 #endif
       
   590     if (segment_within_heart(a, b)) {
       
   591       // Straight line
       
   592       return sqrt((b[0] - a[0]) * (b[0] - a[0]) + (b[1] - a[1]) * (b[1] - a[1]));
       
   593     } else {
       
   594       // Straight line bending around (0, 0)
       
   595       return sqrt(a[0] * a[0] + a[1] * a[1]) + sqrt(b[0] * b[0] + b[1] * b[1]);
       
   596     }
       
   597   }
       
   598 
       
   599   point move_position_toward(point a, double fraction, point b) const 
       
   600   {
       
   601 #ifndef BOOST_NO_STDC_NAMESPACE
       
   602     using std::sqrt;
       
   603 #endif
       
   604 
       
   605     if (segment_within_heart(a, b)) {
       
   606       // Straight line
       
   607       return point(a[0] + (b[0] - a[0]) * fraction,
       
   608                    a[1] + (b[1] - a[1]) * fraction);
       
   609     } else {
       
   610       double distance_to_point_a = sqrt(a[0] * a[0] + a[1] * a[1]);
       
   611       double distance_to_point_b = sqrt(b[0] * b[0] + b[1] * b[1]);
       
   612       double location_of_point = distance_to_point_a / 
       
   613                                    (distance_to_point_a + distance_to_point_b);
       
   614       if (fraction < location_of_point)
       
   615         return point(a[0] * (1 - fraction / location_of_point), 
       
   616                      a[1] * (1 - fraction / location_of_point));
       
   617       else
       
   618         return point(
       
   619           b[0] * ((fraction - location_of_point) / (1 - location_of_point)),
       
   620           b[1] * ((fraction - location_of_point) / (1 - location_of_point)));
       
   621     }
       
   622   }
       
   623 
       
   624  private:
       
   625   shared_ptr<RandomNumberGenerator> gen_ptr;
       
   626   shared_ptr<rand_t> rand;
       
   627 };
       
   628 
       
   629 } // namespace boost
       
   630 
       
   631 #endif // BOOST_GRAPH_GURSOY_ATUN_LAYOUT_HPP