random.hpp

Go to the documentation of this file.

#ifndef _RANDOM_HPP
#define _RANDOM_HPP

#include <iostream>
#include <assert.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "geometry.hpp"

namespace Arak {
  namespace Util {

    class Random {

    protected:

      gsl_rng* r;

    public:

      Random() {
  const gsl_rng_type * T;
  gsl_rng_env_setup();
  T = gsl_rng_default;
  r = gsl_rng_alloc (T);
      }

      Random(const Random& other) : r(gsl_rng_clone(other.r)) { }

      ~Random() {
  gsl_rng_free(r);
      }

      double uniform() { 
  return gsl_rng_uniform(r);
      }

      double uniform(double a, double b) { 
  assert(a < b);
  double x = gsl_rng_uniform(r);
  return x * a + (1 - x) * b;
      }

      unsigned long int uniform(unsigned long int n) { 
  return gsl_rng_uniform_int(r, n);
      }

      bool bernoulli(double p = 0.5) {
  return (uniform() < p);
      }

      unsigned int binomial(double p, unsigned int n) {
  return gsl_ran_binomial(r, p, n);
      }

      inline Geometry::Point uniform(const Geometry::Point& source,
             const Geometry::Point& target) {
  Geometry::Vector a(CGAL::ORIGIN, source);
  Geometry::Vector b(CGAL::ORIGIN, target);
  Geometry::Kernel::FT x(uniform());
  Geometry::Vector u = (a * x) + b - (b * x);
  Geometry::Point p(u.x(), u.y());
#ifdef SANITY_CHECK
  assert(CGAL::collinear_are_ordered_along_line(s.source(), p, s.target()));
#endif
  return p;
      }
      
      inline Geometry::Point uniform(const Geometry::Segment& s) {
  return uniform(s.source(), s.target());
      }

      inline Geometry::Point uniform(const Geometry::Triangle& t) {
  // Generate two uniform random numbers.
  double r1 = uniform();
  double r2 = uniform();
  // Compute the barycentric coordinates of the random point.
  double sqrt_r1 = sqrt(r1);
  double alpha = 1.0 - sqrt_r1;
  double beta = (1.0 - r2) * sqrt_r1;
  double gamma = r2 * sqrt_r1;
  // Translate to Euclidean coordinates.
  return CGAL::ORIGIN 
    + (Geometry::Vector(CGAL::ORIGIN, t[0]) * alpha)
    + (Geometry::Vector(CGAL::ORIGIN, t[1]) * beta)
    + (Geometry::Vector(CGAL::ORIGIN, t[2]) * gamma);
      }

      void write(std::ostream& out) const {
  out << gsl_rng_name(r) 
      << " (with seed " << gsl_rng_default_seed << ")";
      }
    };

    extern Random default_random;

  } // end namespace: Arak::Util
} // end namespace: Arak

inline std::ostream& operator<<(std::ostream& out, 
        const Arak::Util::Random& r) {
  r.write(out);
  return out;
}

#endif

---