boost/graph/topology.hpp

   1 // Copyright 2009 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_TOPOLOGY_HPP
  11 #define BOOST_GRAPH_TOPOLOGY_HPP
  12
  13 #include <boost/config/no_tr1/cmath.hpp>
  14 #include <cmath>
  15 #include <boost/random/uniform_01.hpp>
  16 #include <boost/random/linear_congruential.hpp>
  17 #include <boost/math/constants/constants.hpp> // For root_two
  18 #include <boost/algorithm/minmax.hpp>
  19 #include <boost/config.hpp> // For BOOST_STATIC_CONSTANT
  20 #include <boost/math/special_functions/hypot.hpp>
  21
  22 // Classes and concepts to represent points in a space, with distance and move
  23 // operations (used for Gurson-Atun layout), plus other things like bounding
  24 // boxes used for other layout algorithms.
  25
  26 namespace boost {
  27
  28 /***********************************************************
  29  * Topologies                                              *
  30  ***********************************************************/
  31 template<std::size_t Dims>
  32 class convex_topology
  33 {
  34   public: // For VisualAge C++
  35   struct point
  36   {
  37     BOOST_STATIC_CONSTANT(std::size_t, dimensions = Dims);
  38     point() { }
  39     double& operator[](std::size_t i) {return values[i];}
  40     const double& operator[](std::size_t i) const {return values[i];}
  41
  42   private:
  43     double values[Dims];
  44   };
  45
  46   public: // For VisualAge C++
  47   struct point_difference
  48   {
  49     BOOST_STATIC_CONSTANT(std::size_t, dimensions = Dims);
  50     point_difference() {
  51       for (std::size_t i = 0; i < Dims; ++i) values[i] = 0.;
  52     }
  53     double& operator[](std::size_t i) {return values[i];}
  54     const double& operator[](std::size_t i) const {return values[i];}
  55
  56     friend point_difference operator+(const point_difference& a, const point_difference& b) {
  57       point_difference result;
  58       for (std::size_t i = 0; i < Dims; ++i)
  59         result[i] = a[i] + b[i];
  60       return result;
  61     }
  62
  63     friend point_difference& operator+=(point_difference& a, const point_difference& b) {
  64       for (std::size_t i = 0; i < Dims; ++i)
  65         a[i] += b[i];
  66       return a;
  67     }
  68
  69     friend point_difference operator-(const point_difference& a) {
  70       point_difference result;
  71       for (std::size_t i = 0; i < Dims; ++i)
  72         result[i] = -a[i];
  73       return result;
  74     }
  75
  76     friend point_difference operator-(const point_difference& a, const point_difference& b) {
  77       point_difference result;
  78       for (std::size_t i = 0; i < Dims; ++i)
  79         result[i] = a[i] - b[i];
  80       return result;
  81     }
  82
  83     friend point_difference& operator-=(point_difference& a, const point_difference& b) {
  84       for (std::size_t i = 0; i < Dims; ++i)
  85         a[i] -= b[i];
  86       return a;
  87     }
  88
  89     friend point_difference operator*(const point_difference& a, const point_difference& b) {
  90       point_difference result;
  91       for (std::size_t i = 0; i < Dims; ++i)
  92         result[i] = a[i] * b[i];
  93       return result;
  94     }
  95
  96     friend point_difference operator*(const point_difference& a, double b) {
  97       point_difference result;
  98       for (std::size_t i = 0; i < Dims; ++i)
  99         result[i] = a[i] * b;
 100       return result;
 101     }
 102
 103     friend point_difference operator*(double a, const point_difference& b) {
 104       point_difference result;
 105       for (std::size_t i = 0; i < Dims; ++i)
 106         result[i] = a * b[i];
 107       return result;
 108     }
 109
 110     friend point_difference operator/(const point_difference& a, const point_difference& b) {
 111       point_difference result;
 112       for (std::size_t i = 0; i < Dims; ++i)
 113         result[i] = (b[i] == 0.) ? 0. : a[i] / b[i];
 114       return result;
 115     }
 116
 117     friend double dot(const point_difference& a, const point_difference& b) {
 118       double result = 0;
 119       for (std::size_t i = 0; i < Dims; ++i)
 120         result += a[i] * b[i];
 121       return result;
 122     }
 123
 124   private:
 125     double values[Dims];
 126   };
 127
 128  public:
 129   typedef point point_type;
 130   typedef point_difference point_difference_type;
 131
 132   double distance(point a, point b) const
 133   {
 134     double dist = 0.;
 135     for (std::size_t i = 0; i < Dims; ++i) {
 136       double diff = b[i] - a[i];
 137       dist = boost::math::hypot(dist, diff);
 138     }
 139     // Exact properties of the distance are not important, as long as
 140     // < on what this returns matches real distances; l_2 is used because
 141     // Fruchterman-Reingold also uses this code and it relies on l_2.
 142     return dist;
 143   }
 144
 145   point move_position_toward(point a, double fraction, point b) const
 146   {
 147     point result;
 148     for (std::size_t i = 0; i < Dims; ++i)
 149       result[i] = a[i] + (b[i] - a[i]) * fraction;
 150     return result;
 151   }
 152
 153   point_difference difference(point a, point b) const {
 154     point_difference result;
 155     for (std::size_t i = 0; i < Dims; ++i)
 156       result[i] = a[i] - b[i];
 157     return result;
 158   }
 159
 160   point adjust(point a, point_difference delta) const {
 161     point result;
 162     for (std::size_t i = 0; i < Dims; ++i)
 163       result[i] = a[i] + delta[i];
 164     return result;
 165   }
 166
 167   point pointwise_min(point a, point b) const {
 168     BOOST_USING_STD_MIN();
 169     point result;
 170     for (std::size_t i = 0; i < Dims; ++i)
 171       result[i] = min BOOST_PREVENT_MACRO_SUBSTITUTION (a[i], b[i]);
 172     return result;
 173   }
 174
 175   point pointwise_max(point a, point b) const {
 176     BOOST_USING_STD_MAX();
 177     point result;
 178     for (std::size_t i = 0; i < Dims; ++i)
 179       result[i] = max BOOST_PREVENT_MACRO_SUBSTITUTION (a[i], b[i]);
 180     return result;
 181   }
 182
 183   double norm(point_difference delta) const {
 184     double n = 0.;
 185     for (std::size_t i = 0; i < Dims; ++i)
 186       n = boost::math::hypot(n, delta[i]);
 187     return n;
 188   }
 189
 190   double volume(point_difference delta) const {
 191     double n = 1.;
 192     for (std::size_t i = 0; i < Dims; ++i)
 193       n *= delta[i];
 194     return n;
 195   }
 196
 197 };
 198
 199 template<std::size_t Dims,
 200          typename RandomNumberGenerator = minstd_rand>
 201 class hypercube_topology : public convex_topology<Dims>
 202 {
 203   typedef uniform_01<RandomNumberGenerator, double> rand_t;
 204
 205  public:
 206   typedef typename convex_topology<Dims>::point_type point_type;
 207   typedef typename convex_topology<Dims>::point_difference_type point_difference_type;
 208
 209   explicit hypercube_topology(double scaling = 1.0)
 210     : gen_ptr(new RandomNumberGenerator), rand(new rand_t(*gen_ptr)),
 211       scaling(scaling)
 212   { }
 213
 214   hypercube_topology(RandomNumberGenerator& gen, double scaling = 1.0)
 215     : gen_ptr(), rand(new rand_t(gen)), scaling(scaling) { }
 216
 217   point_type random_point() const
 218   {
 219     point_type p;
 220     for (std::size_t i = 0; i < Dims; ++i)
 221       p[i] = (*rand)() * scaling;
 222     return p;
 223   }
 224
 225   point_type bound(point_type a) const
 226   {
 227     BOOST_USING_STD_MIN();
 228     BOOST_USING_STD_MAX();
 229     point_type p;
 230     for (std::size_t i = 0; i < Dims; ++i)
 231       p[i] = min BOOST_PREVENT_MACRO_SUBSTITUTION (scaling, max BOOST_PREVENT_MACRO_SUBSTITUTION (-scaling, a[i]));
 232     return p;
 233   }
 234
 235   double distance_from_boundary(point_type a) const
 236   {
 237     BOOST_USING_STD_MIN();
 238     BOOST_USING_STD_MAX();
 239 #ifndef BOOST_NO_STDC_NAMESPACE
 240     using std::abs;
 241 #endif
 242     BOOST_STATIC_ASSERT (Dims >= 1);
 243     double dist = abs(scaling - a[0]);
 244     for (std::size_t i = 1; i < Dims; ++i)
 245       dist = min BOOST_PREVENT_MACRO_SUBSTITUTION (dist, abs(scaling - a[i]));
 246     return dist;
 247   }
 248
 249   point_type center() const {
 250     point_type result;
 251     for (std::size_t i = 0; i < Dims; ++i)
 252       result[i] = 0;
 253     return result;
 254   }
 255
 256  private:
 257   shared_ptr<RandomNumberGenerator> gen_ptr;
 258   shared_ptr<rand_t> rand;
 259   double scaling;
 260 };
 261
 262 template<typename RandomNumberGenerator = minstd_rand>
 263 class square_topology : public hypercube_topology<2, RandomNumberGenerator>
 264 {
 265   typedef hypercube_topology<2, RandomNumberGenerator> inherited;
 266
 267  public:
 268   explicit square_topology(double scaling = 1.0) : inherited(scaling) { }
 269
 270   square_topology(RandomNumberGenerator& gen, double scaling = 1.0)
 271     : inherited(gen, scaling) { }
 272 };
 273
 274 template<typename RandomNumberGenerator = minstd_rand>
 275 class rectangle_topology : public convex_topology<2>
 276 {
 277   typedef uniform_01<RandomNumberGenerator, double> rand_t;
 278
 279   public:
 280   rectangle_topology(double left, double top, double right, double bottom)
 281     : gen_ptr(new RandomNumberGenerator), rand(new rand_t(*gen_ptr)),
 282       left(std::min BOOST_PREVENT_MACRO_SUBSTITUTION (left, right)),
 283       top(std::min BOOST_PREVENT_MACRO_SUBSTITUTION (top, bottom)),
 284       right(std::max BOOST_PREVENT_MACRO_SUBSTITUTION (left, right)),
 285       bottom(std::max BOOST_PREVENT_MACRO_SUBSTITUTION (top, bottom)) { }
 286
 287   rectangle_topology(RandomNumberGenerator& gen, double left, double top, double right, double bottom)
 288     : gen_ptr(), rand(new rand_t(gen)),
 289       left(std::min BOOST_PREVENT_MACRO_SUBSTITUTION (left, right)),
 290       top(std::min BOOST_PREVENT_MACRO_SUBSTITUTION (top, bottom)),
 291       right(std::max BOOST_PREVENT_MACRO_SUBSTITUTION (left, right)),
 292       bottom(std::max BOOST_PREVENT_MACRO_SUBSTITUTION (top, bottom)) { }
 293
 294   typedef typename convex_topology<2>::point_type point_type;
 295   typedef typename convex_topology<2>::point_difference_type point_difference_type;
 296
 297   point_type random_point() const
 298   {
 299     point_type p;
 300     p[0] = (*rand)() * (right - left) + left;
 301     p[1] = (*rand)() * (bottom - top) + top;
 302     return p;
 303   }
 304
 305   point_type bound(point_type a) const
 306   {
 307     BOOST_USING_STD_MIN();
 308     BOOST_USING_STD_MAX();
 309     point_type p;
 310     p[0] = min BOOST_PREVENT_MACRO_SUBSTITUTION (right, max BOOST_PREVENT_MACRO_SUBSTITUTION (left, a[0]));
 311     p[1] = min BOOST_PREVENT_MACRO_SUBSTITUTION (bottom, max BOOST_PREVENT_MACRO_SUBSTITUTION (top, a[1]));
 312     return p;
 313   }
 314
 315   double distance_from_boundary(point_type a) const
 316   {
 317     BOOST_USING_STD_MIN();
 318     BOOST_USING_STD_MAX();
 319 #ifndef BOOST_NO_STDC_NAMESPACE
 320     using std::abs;
 321 #endif
 322     double dist = abs(left - a[0]);
 323     dist = min BOOST_PREVENT_MACRO_SUBSTITUTION (dist, abs(right - a[0]));
 324     dist = min BOOST_PREVENT_MACRO_SUBSTITUTION (dist, abs(top - a[1]));
 325     dist = min BOOST_PREVENT_MACRO_SUBSTITUTION (dist, abs(bottom - a[1]));
 326     return dist;
 327   }
 328
 329   point_type center() const {
 330     point_type result;
 331     result[0] = (left + right) / 2.;
 332     result[1] = (top + bottom) / 2.;
 333     return result;
 334   }
 335
 336  private:
 337   shared_ptr<RandomNumberGenerator> gen_ptr;
 338   shared_ptr<rand_t> rand;
 339   double left, top, right, bottom;
 340 };
 341
 342 template<typename RandomNumberGenerator = minstd_rand>
 343 class cube_topology : public hypercube_topology<3, RandomNumberGenerator>
 344 {
 345   typedef hypercube_topology<3, RandomNumberGenerator> inherited;
 346
 347  public:
 348   explicit cube_topology(double scaling = 1.0) : inherited(scaling) { }
 349
 350   cube_topology(RandomNumberGenerator& gen, double scaling = 1.0)
 351     : inherited(gen, scaling) { }
 352 };
 353
 354 template<std::size_t Dims,
 355          typename RandomNumberGenerator = minstd_rand>
 356 class ball_topology : public convex_topology<Dims>
 357 {
 358   typedef uniform_01<RandomNumberGenerator, double> rand_t;
 359
 360  public:
 361   typedef typename convex_topology<Dims>::point_type point_type;
 362
 363   explicit ball_topology(double radius = 1.0)
 364     : gen_ptr(new RandomNumberGenerator), rand(new rand_t(*gen_ptr)),
 365       radius(radius)
 366   { }
 367
 368   ball_topology(RandomNumberGenerator& gen, double radius = 1.0)
 369     : gen_ptr(), rand(new rand_t(gen)), radius(radius) { }
 370
 371   point_type random_point() const
 372   {
 373     point_type p;
 374     double dist_sum;
 375     do {
 376       dist_sum = 0.0;
 377       for (std::size_t i = 0; i < Dims; ++i) {
 378         double x = (*rand)() * 2*radius - radius;
 379         p[i] = x;
 380         dist_sum += x * x;
 381       }
 382     } while (dist_sum > radius*radius);
 383     return p;
 384   }
 385
 386   point_type bound(point_type a) const
 387   {
 388     BOOST_USING_STD_MIN();
 389     BOOST_USING_STD_MAX();
 390     double r = 0.;
 391     for (std::size_t i = 0; i < Dims; ++i)
 392       r = boost::math::hypot(r, a[i]);
 393     if (r <= radius) return a;
 394     double scaling_factor = radius / r;
 395     point_type p;
 396     for (std::size_t i = 0; i < Dims; ++i)
 397       p[i] = a[i] * scaling_factor;
 398     return p;
 399   }
 400
 401   double distance_from_boundary(point_type a) const
 402   {
 403     double r = 0.;
 404     for (std::size_t i = 0; i < Dims; ++i)
 405       r = boost::math::hypot(r, a[i]);
 406     return radius - r;
 407   }
 408
 409   point_type center() const {
 410     point_type result;
 411     for (std::size_t i = 0; i < Dims; ++i)
 412       result[i] = 0;
 413     return result;
 414   }
 415
 416  private:
 417   shared_ptr<RandomNumberGenerator> gen_ptr;
 418   shared_ptr<rand_t> rand;
 419   double radius;
 420 };
 421
 422 template<typename RandomNumberGenerator = minstd_rand>
 423 class circle_topology : public ball_topology<2, RandomNumberGenerator>
 424 {
 425   typedef ball_topology<2, RandomNumberGenerator> inherited;
 426
 427  public:
 428   explicit circle_topology(double radius = 1.0) : inherited(radius) { }
 429
 430   circle_topology(RandomNumberGenerator& gen, double radius = 1.0)
 431     : inherited(gen, radius) { }
 432 };
 433
 434 template<typename RandomNumberGenerator = minstd_rand>
 435 class sphere_topology : public ball_topology<3, RandomNumberGenerator>
 436 {
 437   typedef ball_topology<3, RandomNumberGenerator> inherited;
 438
 439  public:
 440   explicit sphere_topology(double radius = 1.0) : inherited(radius) { }
 441
 442   sphere_topology(RandomNumberGenerator& gen, double radius = 1.0)
 443     : inherited(gen, radius) { }
 444 };
 445
 446 template<typename RandomNumberGenerator = minstd_rand>
 447 class heart_topology
 448 {
 449   // Heart is defined as the union of three shapes:
 450   // Square w/ corners (+-1000, -1000), (0, 0), (0, -2000)
 451   // Circle centered at (-500, -500) radius 500*sqrt(2)
 452   // Circle centered at (500, -500) radius 500*sqrt(2)
 453   // Bounding box (-1000, -2000) - (1000, 500*(sqrt(2) - 1))
 454
 455   struct point
 456   {
 457     point() { values[0] = 0.0; values[1] = 0.0; }
 458     point(double x, double y) { values[0] = x; values[1] = y; }
 459
 460     double& operator[](std::size_t i)       { return values[i]; }
 461     double  operator[](std::size_t i) const { return values[i]; }
 462
 463   private:
 464     double values[2];
 465   };
 466
 467   bool in_heart(point p) const
 468   {
 469 #ifndef BOOST_NO_STDC_NAMESPACE
 470     using std::abs;
 471 #endif
 472
 473     if (p[1] < abs(p[0]) - 2000) return false; // Bottom
 474     if (p[1] <= -1000) return true; // Diagonal of square
 475     if (boost::math::hypot(p[0] - -500, p[1] - -500) <= 500. * boost::math::constants::root_two<double>())
 476       return true; // Left circle
 477     if (boost::math::hypot(p[0] - 500, p[1] - -500) <= 500. * boost::math::constants::root_two<double>())
 478       return true; // Right circle
 479     return false;
 480   }
 481
 482   bool segment_within_heart(point p1, point p2) const
 483   {
 484     // Assumes that p1 and p2 are within the heart
 485     if ((p1[0] < 0) == (p2[0] < 0)) return true; // Same side of symmetry line
 486     if (p1[0] == p2[0]) return true; // Vertical
 487     double slope = (p2[1] - p1[1]) / (p2[0] - p1[0]);
 488     double intercept = p1[1] - p1[0] * slope;
 489     if (intercept > 0) return false; // Crosses between circles
 490     return true;
 491   }
 492
 493   typedef uniform_01<RandomNumberGenerator, double> rand_t;
 494
 495  public:
 496   typedef point point_type;
 497
 498   heart_topology()
 499     : gen_ptr(new RandomNumberGenerator), rand(new rand_t(*gen_ptr)) { }
 500
 501   heart_topology(RandomNumberGenerator& gen)
 502     : gen_ptr(), rand(new rand_t(gen)) { }
 503
 504   point random_point() const
 505   {
 506     point result;
 507     do {
 508       result[0] = (*rand)() * (1000 + 1000 * boost::math::constants::root_two<double>()) - (500 + 500 * boost::math::constants::root_two<double>());
 509       result[1] = (*rand)() * (2000 + 500 * (boost::math::constants::root_two<double>() - 1)) - 2000;
 510     } while (!in_heart(result));
 511     return result;
 512   }
 513
 514   // Not going to provide clipping to bounding region or distance from boundary
 515
 516   double distance(point a, point b) const
 517   {
 518     if (segment_within_heart(a, b)) {
 519       // Straight line
 520       return boost::math::hypot(b[0] - a[0], b[1] - a[1]);
 521     } else {
 522       // Straight line bending around (0, 0)
 523       return boost::math::hypot(a[0], a[1]) + boost::math::hypot(b[0], b[1]);
 524     }
 525   }
 526
 527   point move_position_toward(point a, double fraction, point b) const
 528   {
 529     if (segment_within_heart(a, b)) {
 530       // Straight line
 531       return point(a[0] + (b[0] - a[0]) * fraction,
 532                    a[1] + (b[1] - a[1]) * fraction);
 533     } else {
 534       double distance_to_point_a = boost::math::hypot(a[0], a[1]);
 535       double distance_to_point_b = boost::math::hypot(b[0], b[1]);
 536       double location_of_point = distance_to_point_a /
 537                                    (distance_to_point_a + distance_to_point_b);
 538       if (fraction < location_of_point)
 539         return point(a[0] * (1 - fraction / location_of_point),
 540                      a[1] * (1 - fraction / location_of_point));
 541       else
 542         return point(
 543           b[0] * ((fraction - location_of_point) / (1 - location_of_point)),
 544           b[1] * ((fraction - location_of_point) / (1 - location_of_point)));
 545     }
 546   }
 547
 548  private:
 549   shared_ptr<RandomNumberGenerator> gen_ptr;
 550   shared_ptr<rand_t> rand;
 551 };
 552
 553 } // namespace boost
 554
 555 #endif // BOOST_GRAPH_TOPOLOGY_HPP