Chris@16: // Copyright 2009 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: Jeremiah Willcock
Chris@16: //           Douglas Gregor
Chris@16: //           Andrew Lumsdaine
Chris@16: #ifndef BOOST_GRAPH_TOPOLOGY_HPP
Chris@16: #define BOOST_GRAPH_TOPOLOGY_HPP
Chris@16: 
Chris@16: #include <boost/config/no_tr1/cmath.hpp>
Chris@16: #include <cmath>
Chris@16: #include <boost/random/uniform_01.hpp>
Chris@16: #include <boost/random/linear_congruential.hpp>
Chris@16: #include <boost/math/constants/constants.hpp> // For root_two
Chris@16: #include <boost/algorithm/minmax.hpp>
Chris@16: #include <boost/config.hpp> // For BOOST_STATIC_CONSTANT
Chris@16: #include <boost/math/special_functions/hypot.hpp>
Chris@16: 
Chris@16: // Classes and concepts to represent points in a space, with distance and move
Chris@16: // operations (used for Gurson-Atun layout), plus other things like bounding
Chris@16: // boxes used for other layout algorithms.
Chris@16: 
Chris@16: namespace boost {
Chris@16: 
Chris@16: /***********************************************************
Chris@16:  * Topologies                                              *
Chris@16:  ***********************************************************/
Chris@16: template<std::size_t Dims>
Chris@16: class convex_topology 
Chris@16: {
Chris@16:   public: // For VisualAge C++
Chris@16:   struct point 
Chris@16:   {
Chris@16:     BOOST_STATIC_CONSTANT(std::size_t, dimensions = Dims);
Chris@16:     point() { }
Chris@16:     double& operator[](std::size_t i) {return values[i];}
Chris@16:     const double& operator[](std::size_t i) const {return values[i];}
Chris@16: 
Chris@16:   private:
Chris@16:     double values[Dims];
Chris@16:   };
Chris@16: 
Chris@16:   public: // For VisualAge C++
Chris@16:   struct point_difference
Chris@16:   {
Chris@16:     BOOST_STATIC_CONSTANT(std::size_t, dimensions = Dims);
Chris@16:     point_difference() {
Chris@16:       for (std::size_t i = 0; i < Dims; ++i) values[i] = 0.;
Chris@16:     }
Chris@16:     double& operator[](std::size_t i) {return values[i];}
Chris@16:     const double& operator[](std::size_t i) const {return values[i];}
Chris@16: 
Chris@16:     friend point_difference operator+(const point_difference& a, const point_difference& b) {
Chris@16:       point_difference result;
Chris@16:       for (std::size_t i = 0; i < Dims; ++i)
Chris@16:         result[i] = a[i] + b[i];
Chris@16:       return result;
Chris@16:     }
Chris@16: 
Chris@16:     friend point_difference& operator+=(point_difference& a, const point_difference& b) {
Chris@16:       for (std::size_t i = 0; i < Dims; ++i)
Chris@16:         a[i] += b[i];
Chris@16:       return a;
Chris@16:     }
Chris@16: 
Chris@16:     friend point_difference operator-(const point_difference& a) {
Chris@16:       point_difference result;
Chris@16:       for (std::size_t i = 0; i < Dims; ++i)
Chris@16:         result[i] = -a[i];
Chris@16:       return result;
Chris@16:     }
Chris@16: 
Chris@16:     friend point_difference operator-(const point_difference& a, const point_difference& b) {
Chris@16:       point_difference result;
Chris@16:       for (std::size_t i = 0; i < Dims; ++i)
Chris@16:         result[i] = a[i] - b[i];
Chris@16:       return result;
Chris@16:     }
Chris@16: 
Chris@16:     friend point_difference& operator-=(point_difference& a, const point_difference& b) {
Chris@16:       for (std::size_t i = 0; i < Dims; ++i)
Chris@16:         a[i] -= b[i];
Chris@16:       return a;
Chris@16:     }
Chris@16: 
Chris@16:     friend point_difference operator*(const point_difference& a, const point_difference& b) {
Chris@16:       point_difference result;
Chris@16:       for (std::size_t i = 0; i < Dims; ++i)
Chris@16:         result[i] = a[i] * b[i];
Chris@16:       return result;
Chris@16:     }
Chris@16: 
Chris@16:     friend point_difference operator*(const point_difference& a, double b) {
Chris@16:       point_difference result;
Chris@16:       for (std::size_t i = 0; i < Dims; ++i)
Chris@16:         result[i] = a[i] * b;
Chris@16:       return result;
Chris@16:     }
Chris@16: 
Chris@16:     friend point_difference operator*(double a, const point_difference& b) {
Chris@16:       point_difference result;
Chris@16:       for (std::size_t i = 0; i < Dims; ++i)
Chris@16:         result[i] = a * b[i];
Chris@16:       return result;
Chris@16:     }
Chris@16: 
Chris@16:     friend point_difference operator/(const point_difference& a, const point_difference& b) {
Chris@16:       point_difference result;
Chris@16:       for (std::size_t i = 0; i < Dims; ++i)
Chris@16:         result[i] = (b[i] == 0.) ? 0. : a[i] / b[i];
Chris@16:       return result;
Chris@16:     }
Chris@16: 
Chris@16:     friend double dot(const point_difference& a, const point_difference& b) {
Chris@16:       double result = 0;
Chris@16:       for (std::size_t i = 0; i < Dims; ++i)
Chris@16:         result += a[i] * b[i];
Chris@16:       return result;
Chris@16:     }
Chris@16: 
Chris@16:   private:
Chris@16:     double values[Dims];
Chris@16:   };
Chris@16: 
Chris@16:  public:
Chris@16:   typedef point point_type;
Chris@16:   typedef point_difference point_difference_type;
Chris@16: 
Chris@16:   double distance(point a, point b) const 
Chris@16:   {
Chris@16:     double dist = 0.;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i) {
Chris@16:       double diff = b[i] - a[i];
Chris@16:       dist = boost::math::hypot(dist, diff);
Chris@16:     }
Chris@16:     // Exact properties of the distance are not important, as long as
Chris@16:     // < on what this returns matches real distances; l_2 is used because
Chris@16:     // Fruchterman-Reingold also uses this code and it relies on l_2.
Chris@16:     return dist;
Chris@16:   }
Chris@16: 
Chris@16:   point move_position_toward(point a, double fraction, point b) const 
Chris@16:   {
Chris@16:     point result;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       result[i] = a[i] + (b[i] - a[i]) * fraction;
Chris@16:     return result;
Chris@16:   }
Chris@16: 
Chris@16:   point_difference difference(point a, point b) const {
Chris@16:     point_difference result;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       result[i] = a[i] - b[i];
Chris@16:     return result;
Chris@16:   }
Chris@16: 
Chris@16:   point adjust(point a, point_difference delta) const {
Chris@16:     point result;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       result[i] = a[i] + delta[i];
Chris@16:     return result;
Chris@16:   }
Chris@16: 
Chris@16:   point pointwise_min(point a, point b) const {
Chris@16:     BOOST_USING_STD_MIN();
Chris@16:     point result;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       result[i] = min BOOST_PREVENT_MACRO_SUBSTITUTION (a[i], b[i]);
Chris@16:     return result;
Chris@16:   }
Chris@16: 
Chris@16:   point pointwise_max(point a, point b) const {
Chris@16:     BOOST_USING_STD_MAX();
Chris@16:     point result;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       result[i] = max BOOST_PREVENT_MACRO_SUBSTITUTION (a[i], b[i]);
Chris@16:     return result;
Chris@16:   }
Chris@16: 
Chris@16:   double norm(point_difference delta) const {
Chris@16:     double n = 0.;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       n = boost::math::hypot(n, delta[i]);
Chris@16:     return n;
Chris@16:   }
Chris@16: 
Chris@16:   double volume(point_difference delta) const {
Chris@16:     double n = 1.;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       n *= delta[i];
Chris@16:     return n;
Chris@16:   }
Chris@16: 
Chris@16: };
Chris@16: 
Chris@16: template<std::size_t Dims,
Chris@16:          typename RandomNumberGenerator = minstd_rand>
Chris@16: class hypercube_topology : public convex_topology<Dims>
Chris@16: {
Chris@16:   typedef uniform_01<RandomNumberGenerator, double> rand_t;
Chris@16: 
Chris@16:  public:
Chris@16:   typedef typename convex_topology<Dims>::point_type point_type;
Chris@16:   typedef typename convex_topology<Dims>::point_difference_type point_difference_type;
Chris@16: 
Chris@16:   explicit hypercube_topology(double scaling = 1.0) 
Chris@16:     : gen_ptr(new RandomNumberGenerator), rand(new rand_t(*gen_ptr)), 
Chris@16:       scaling(scaling) 
Chris@16:   { }
Chris@16: 
Chris@16:   hypercube_topology(RandomNumberGenerator& gen, double scaling = 1.0) 
Chris@16:     : gen_ptr(), rand(new rand_t(gen)), scaling(scaling) { }
Chris@16:                      
Chris@16:   point_type random_point() const 
Chris@16:   {
Chris@16:     point_type p;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       p[i] = (*rand)() * scaling;
Chris@16:     return p;
Chris@16:   }
Chris@16: 
Chris@16:   point_type bound(point_type a) const
Chris@16:   {
Chris@16:     BOOST_USING_STD_MIN();
Chris@16:     BOOST_USING_STD_MAX();
Chris@16:     point_type p;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       p[i] = min BOOST_PREVENT_MACRO_SUBSTITUTION (scaling, max BOOST_PREVENT_MACRO_SUBSTITUTION (-scaling, a[i]));
Chris@16:     return p;
Chris@16:   }
Chris@16: 
Chris@16:   double distance_from_boundary(point_type a) const
Chris@16:   {
Chris@16:     BOOST_USING_STD_MIN();
Chris@16:     BOOST_USING_STD_MAX();
Chris@16: #ifndef BOOST_NO_STDC_NAMESPACE
Chris@16:     using std::abs;
Chris@16: #endif
Chris@16:     BOOST_STATIC_ASSERT (Dims >= 1);
Chris@16:     double dist = abs(scaling - a[0]);
Chris@16:     for (std::size_t i = 1; i < Dims; ++i)
Chris@16:       dist = min BOOST_PREVENT_MACRO_SUBSTITUTION (dist, abs(scaling - a[i]));
Chris@16:     return dist;
Chris@16:   }
Chris@16: 
Chris@16:   point_type center() const {
Chris@16:     point_type result;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       result[i] = scaling * .5;
Chris@16:     return result;
Chris@16:   }
Chris@16: 
Chris@16:   point_type origin() const {
Chris@16:     point_type result;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       result[i] = 0;
Chris@16:     return result;
Chris@16:   }
Chris@16: 
Chris@16:   point_difference_type extent() const {
Chris@16:     point_difference_type result;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       result[i] = scaling;
Chris@16:     return result;
Chris@16:   }
Chris@16: 
Chris@16:  private:
Chris@16:   shared_ptr<RandomNumberGenerator> gen_ptr;
Chris@16:   shared_ptr<rand_t> rand;
Chris@16:   double scaling;
Chris@16: };
Chris@16: 
Chris@16: template<typename RandomNumberGenerator = minstd_rand>
Chris@16: class square_topology : public hypercube_topology<2, RandomNumberGenerator>
Chris@16: {
Chris@16:   typedef hypercube_topology<2, RandomNumberGenerator> inherited;
Chris@16: 
Chris@16:  public:
Chris@16:   explicit square_topology(double scaling = 1.0) : inherited(scaling) { }
Chris@16:   
Chris@16:   square_topology(RandomNumberGenerator& gen, double scaling = 1.0) 
Chris@16:     : inherited(gen, scaling) { }
Chris@16: };
Chris@16: 
Chris@16: template<typename RandomNumberGenerator = minstd_rand>
Chris@16: class rectangle_topology : public convex_topology<2>
Chris@16: {
Chris@16:   typedef uniform_01<RandomNumberGenerator, double> rand_t;
Chris@16: 
Chris@16:   public:
Chris@16:   rectangle_topology(double left, double top, double right, double bottom)
Chris@16:     : gen_ptr(new RandomNumberGenerator), rand(new rand_t(*gen_ptr)),
Chris@16:       left(std::min BOOST_PREVENT_MACRO_SUBSTITUTION (left, right)),
Chris@16:       top(std::min BOOST_PREVENT_MACRO_SUBSTITUTION (top, bottom)),
Chris@16:       right(std::max BOOST_PREVENT_MACRO_SUBSTITUTION (left, right)),
Chris@16:       bottom(std::max BOOST_PREVENT_MACRO_SUBSTITUTION (top, bottom)) { }
Chris@16: 
Chris@16:   rectangle_topology(RandomNumberGenerator& gen, double left, double top, double right, double bottom)
Chris@16:     : gen_ptr(), rand(new rand_t(gen)),
Chris@16:       left(std::min BOOST_PREVENT_MACRO_SUBSTITUTION (left, right)),
Chris@16:       top(std::min BOOST_PREVENT_MACRO_SUBSTITUTION (top, bottom)),
Chris@16:       right(std::max BOOST_PREVENT_MACRO_SUBSTITUTION (left, right)),
Chris@16:       bottom(std::max BOOST_PREVENT_MACRO_SUBSTITUTION (top, bottom)) { }
Chris@16: 
Chris@16:   typedef typename convex_topology<2>::point_type point_type;
Chris@16:   typedef typename convex_topology<2>::point_difference_type point_difference_type;
Chris@16: 
Chris@16:   point_type random_point() const 
Chris@16:   {
Chris@16:     point_type p;
Chris@16:     p[0] = (*rand)() * (right - left) + left;
Chris@16:     p[1] = (*rand)() * (bottom - top) + top;
Chris@16:     return p;
Chris@16:   }
Chris@16: 
Chris@16:   point_type bound(point_type a) const
Chris@16:   {
Chris@16:     BOOST_USING_STD_MIN();
Chris@16:     BOOST_USING_STD_MAX();
Chris@16:     point_type p;
Chris@16:     p[0] = min BOOST_PREVENT_MACRO_SUBSTITUTION (right, max BOOST_PREVENT_MACRO_SUBSTITUTION (left, a[0]));
Chris@16:     p[1] = min BOOST_PREVENT_MACRO_SUBSTITUTION (bottom, max BOOST_PREVENT_MACRO_SUBSTITUTION (top, a[1]));
Chris@16:     return p;
Chris@16:   }
Chris@16: 
Chris@16:   double distance_from_boundary(point_type a) const
Chris@16:   {
Chris@16:     BOOST_USING_STD_MIN();
Chris@16:     BOOST_USING_STD_MAX();
Chris@16: #ifndef BOOST_NO_STDC_NAMESPACE
Chris@16:     using std::abs;
Chris@16: #endif
Chris@16:     double dist = abs(left - a[0]);
Chris@16:     dist = min BOOST_PREVENT_MACRO_SUBSTITUTION (dist, abs(right - a[0]));
Chris@16:     dist = min BOOST_PREVENT_MACRO_SUBSTITUTION (dist, abs(top - a[1]));
Chris@16:     dist = min BOOST_PREVENT_MACRO_SUBSTITUTION (dist, abs(bottom - a[1]));
Chris@16:     return dist;
Chris@16:   }
Chris@16: 
Chris@16:   point_type center() const {
Chris@16:     point_type result;
Chris@16:     result[0] = (left + right) / 2.;
Chris@16:     result[1] = (top + bottom) / 2.;
Chris@16:     return result;
Chris@16:   }
Chris@16: 
Chris@16:   point_type origin() const {
Chris@16:     point_type result;
Chris@16:     result[0] = left;
Chris@16:     result[1] = top;
Chris@16:     return result;
Chris@16:   }
Chris@16: 
Chris@16:   point_difference_type extent() const {
Chris@16:     point_difference_type result;
Chris@16:     result[0] = right - left;
Chris@16:     result[1] = bottom - top;
Chris@16:     return result;
Chris@16:   }
Chris@16: 
Chris@16:  private:
Chris@16:   shared_ptr<RandomNumberGenerator> gen_ptr;
Chris@16:   shared_ptr<rand_t> rand;
Chris@16:   double left, top, right, bottom;
Chris@16: };
Chris@16: 
Chris@16: template<typename RandomNumberGenerator = minstd_rand>
Chris@16: class cube_topology : public hypercube_topology<3, RandomNumberGenerator>
Chris@16: {
Chris@16:   typedef hypercube_topology<3, RandomNumberGenerator> inherited;
Chris@16: 
Chris@16:  public:
Chris@16:   explicit cube_topology(double scaling = 1.0) : inherited(scaling) { }
Chris@16:   
Chris@16:   cube_topology(RandomNumberGenerator& gen, double scaling = 1.0) 
Chris@16:     : inherited(gen, scaling) { }
Chris@16: };
Chris@16: 
Chris@16: template<std::size_t Dims,
Chris@16:          typename RandomNumberGenerator = minstd_rand>
Chris@16: class ball_topology : public convex_topology<Dims>
Chris@16: {
Chris@16:   typedef uniform_01<RandomNumberGenerator, double> rand_t;
Chris@16: 
Chris@16:  public:
Chris@16:   typedef typename convex_topology<Dims>::point_type point_type;
Chris@16:   typedef typename convex_topology<Dims>::point_difference_type point_difference_type;
Chris@16: 
Chris@16:   explicit ball_topology(double radius = 1.0) 
Chris@16:     : gen_ptr(new RandomNumberGenerator), rand(new rand_t(*gen_ptr)), 
Chris@16:       radius(radius) 
Chris@16:   { }
Chris@16: 
Chris@16:   ball_topology(RandomNumberGenerator& gen, double radius = 1.0) 
Chris@16:     : gen_ptr(), rand(new rand_t(gen)), radius(radius) { }
Chris@16:                      
Chris@16:   point_type random_point() const 
Chris@16:   {
Chris@16:     point_type p;
Chris@16:     double dist_sum;
Chris@16:     do {
Chris@16:       dist_sum = 0.0;
Chris@16:       for (std::size_t i = 0; i < Dims; ++i) {
Chris@16:         double x = (*rand)() * 2*radius - radius;
Chris@16:         p[i] = x;
Chris@16:         dist_sum += x * x;
Chris@16:       }
Chris@16:     } while (dist_sum > radius*radius);
Chris@16:     return p;
Chris@16:   }
Chris@16: 
Chris@16:   point_type bound(point_type a) const
Chris@16:   {
Chris@16:     BOOST_USING_STD_MIN();
Chris@16:     BOOST_USING_STD_MAX();
Chris@16:     double r = 0.;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       r = boost::math::hypot(r, a[i]);
Chris@16:     if (r <= radius) return a;
Chris@16:     double scaling_factor = radius / r;
Chris@16:     point_type p;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       p[i] = a[i] * scaling_factor;
Chris@16:     return p;
Chris@16:   }
Chris@16: 
Chris@16:   double distance_from_boundary(point_type a) const
Chris@16:   {
Chris@16:     double r = 0.;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       r = boost::math::hypot(r, a[i]);
Chris@16:     return radius - r;
Chris@16:   }
Chris@16: 
Chris@16:   point_type center() const {
Chris@16:     point_type result;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       result[i] = 0;
Chris@16:     return result;
Chris@16:   }
Chris@16: 
Chris@16:   point_type origin() const {
Chris@16:     point_type result;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       result[i] = -radius;
Chris@16:     return result;
Chris@16:   }
Chris@16: 
Chris@16:   point_difference_type extent() const {
Chris@16:     point_difference_type result;
Chris@16:     for (std::size_t i = 0; i < Dims; ++i)
Chris@16:       result[i] = 2. * radius;
Chris@16:     return result;
Chris@16:   }
Chris@16: 
Chris@16:  private:
Chris@16:   shared_ptr<RandomNumberGenerator> gen_ptr;
Chris@16:   shared_ptr<rand_t> rand;
Chris@16:   double radius;
Chris@16: };
Chris@16: 
Chris@16: template<typename RandomNumberGenerator = minstd_rand>
Chris@16: class circle_topology : public ball_topology<2, RandomNumberGenerator>
Chris@16: {
Chris@16:   typedef ball_topology<2, RandomNumberGenerator> inherited;
Chris@16: 
Chris@16:  public:
Chris@16:   explicit circle_topology(double radius = 1.0) : inherited(radius) { }
Chris@16:   
Chris@16:   circle_topology(RandomNumberGenerator& gen, double radius = 1.0) 
Chris@16:     : inherited(gen, radius) { }
Chris@16: };
Chris@16: 
Chris@16: template<typename RandomNumberGenerator = minstd_rand>
Chris@16: class sphere_topology : public ball_topology<3, RandomNumberGenerator>
Chris@16: {
Chris@16:   typedef ball_topology<3, RandomNumberGenerator> inherited;
Chris@16: 
Chris@16:  public:
Chris@16:   explicit sphere_topology(double radius = 1.0) : inherited(radius) { }
Chris@16:   
Chris@16:   sphere_topology(RandomNumberGenerator& gen, double radius = 1.0) 
Chris@16:     : inherited(gen, radius) { }
Chris@16: };
Chris@16: 
Chris@16: template<typename RandomNumberGenerator = minstd_rand>
Chris@16: class heart_topology 
Chris@16: {
Chris@16:   // Heart is defined as the union of three shapes:
Chris@16:   // Square w/ corners (+-1000, -1000), (0, 0), (0, -2000)
Chris@16:   // Circle centered at (-500, -500) radius 500*sqrt(2)
Chris@16:   // Circle centered at (500, -500) radius 500*sqrt(2)
Chris@16:   // Bounding box (-1000, -2000) - (1000, 500*(sqrt(2) - 1))
Chris@16: 
Chris@16:   struct point 
Chris@16:   {
Chris@16:     point() { values[0] = 0.0; values[1] = 0.0; }
Chris@16:     point(double x, double y) { values[0] = x; values[1] = y; }
Chris@16: 
Chris@16:     double& operator[](std::size_t i)       { return values[i]; }
Chris@16:     double  operator[](std::size_t i) const { return values[i]; }
Chris@16: 
Chris@16:   private:
Chris@16:     double values[2];
Chris@16:   };
Chris@16: 
Chris@16:   bool in_heart(point p) const 
Chris@16:   {
Chris@16: #ifndef BOOST_NO_STDC_NAMESPACE
Chris@16:     using std::abs;
Chris@16: #endif
Chris@16: 
Chris@16:     if (p[1] < abs(p[0]) - 2000) return false; // Bottom
Chris@16:     if (p[1] <= -1000) return true; // Diagonal of square
Chris@16:     if (boost::math::hypot(p[0] - -500, p[1] - -500) <= 500. * boost::math::constants::root_two<double>())
Chris@16:       return true; // Left circle
Chris@16:     if (boost::math::hypot(p[0] - 500, p[1] - -500) <= 500. * boost::math::constants::root_two<double>())
Chris@16:       return true; // Right circle
Chris@16:     return false;
Chris@16:   }
Chris@16: 
Chris@16:   bool segment_within_heart(point p1, point p2) const 
Chris@16:   {
Chris@16:     // Assumes that p1 and p2 are within the heart
Chris@16:     if ((p1[0] < 0) == (p2[0] < 0)) return true; // Same side of symmetry line
Chris@16:     if (p1[0] == p2[0]) return true; // Vertical
Chris@16:     double slope = (p2[1] - p1[1]) / (p2[0] - p1[0]);
Chris@16:     double intercept = p1[1] - p1[0] * slope;
Chris@16:     if (intercept > 0) return false; // Crosses between circles
Chris@16:     return true;
Chris@16:   }
Chris@16: 
Chris@16:   typedef uniform_01<RandomNumberGenerator, double> rand_t;
Chris@16: 
Chris@16:  public:
Chris@16:   typedef point point_type;
Chris@16: 
Chris@16:   heart_topology() 
Chris@16:     : gen_ptr(new RandomNumberGenerator), rand(new rand_t(*gen_ptr)) { }
Chris@16: 
Chris@16:   heart_topology(RandomNumberGenerator& gen) 
Chris@16:     : gen_ptr(), rand(new rand_t(gen)) { }
Chris@16: 
Chris@16:   point random_point() const 
Chris@16:   {
Chris@16:     point result;
Chris@16:     do {
Chris@16:       result[0] = (*rand)() * (1000 + 1000 * boost::math::constants::root_two<double>()) - (500 + 500 * boost::math::constants::root_two<double>());
Chris@16:       result[1] = (*rand)() * (2000 + 500 * (boost::math::constants::root_two<double>() - 1)) - 2000;
Chris@16:     } while (!in_heart(result));
Chris@16:     return result;
Chris@16:   }
Chris@16: 
Chris@16:   // Not going to provide clipping to bounding region or distance from boundary
Chris@16: 
Chris@16:   double distance(point a, point b) const 
Chris@16:   {
Chris@16:     if (segment_within_heart(a, b)) {
Chris@16:       // Straight line
Chris@16:       return boost::math::hypot(b[0] - a[0], b[1] - a[1]);
Chris@16:     } else {
Chris@16:       // Straight line bending around (0, 0)
Chris@16:       return boost::math::hypot(a[0], a[1]) + boost::math::hypot(b[0], b[1]);
Chris@16:     }
Chris@16:   }
Chris@16: 
Chris@16:   point move_position_toward(point a, double fraction, point b) const 
Chris@16:   {
Chris@16:     if (segment_within_heart(a, b)) {
Chris@16:       // Straight line
Chris@16:       return point(a[0] + (b[0] - a[0]) * fraction,
Chris@16:                    a[1] + (b[1] - a[1]) * fraction);
Chris@16:     } else {
Chris@16:       double distance_to_point_a = boost::math::hypot(a[0], a[1]);
Chris@16:       double distance_to_point_b = boost::math::hypot(b[0], b[1]);
Chris@16:       double location_of_point = distance_to_point_a / 
Chris@16:                                    (distance_to_point_a + distance_to_point_b);
Chris@16:       if (fraction < location_of_point)
Chris@16:         return point(a[0] * (1 - fraction / location_of_point), 
Chris@16:                      a[1] * (1 - fraction / location_of_point));
Chris@16:       else
Chris@16:         return point(
Chris@16:           b[0] * ((fraction - location_of_point) / (1 - location_of_point)),
Chris@16:           b[1] * ((fraction - location_of_point) / (1 - location_of_point)));
Chris@16:     }
Chris@16:   }
Chris@16: 
Chris@16:  private:
Chris@16:   shared_ptr<RandomNumberGenerator> gen_ptr;
Chris@16:   shared_ptr<rand_t> rand;
Chris@16: };
Chris@16: 
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #endif // BOOST_GRAPH_TOPOLOGY_HPP