random.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2014-2018
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */

#ifndef SRC_INCLUDE_RANDOM_H_
#define SRC_INCLUDE_RANDOM_H_

#include <cassert>
#include <limits>
#include <random>
#include <utility>
#include <vector>

namespace smash {

namespace random {

using Engine = std::mt19937_64;

extern /*thread_local (see #3075)*/ Engine engine;

template <typename T>
class uniform_dist {
 public:
  uniform_dist(T min, T max) : distribution(min, max) {}
  T operator()() { return distribution(engine); }

 private:
  std::uniform_real_distribution<T> distribution;
};

template <typename T>
void set_seed(T &&seed) {
  static_assert(std::is_same<Engine::result_type, uint64_t>::value,
                "experiment.cc needs the seed to be 64 bits");
  engine.seed(std::forward<T>(seed));
}

inline Engine::result_type advance() { return engine(); }

template <typename T>
T uniform(T min, T max) {
  return std::uniform_real_distribution<T>(min, max)(engine);
}

template <typename T>
T uniform_int(T min, T max) {
  return std::uniform_int_distribution<T>(min, max)(engine);
}

template <typename T = double>
T canonical() {
  return std::generate_canonical<T, std::numeric_limits<double>::digits>(
      engine);
}

template <typename T = double>
T canonical_nonzero() {
  // use 'nextafter' to generate a value that is guaranteed to be larger than 0
  return std::nextafter(
      std::generate_canonical<T, std::numeric_limits<double>::digits>(engine),
      T(1));
}

template <typename T>
uniform_dist<T> make_uniform_distribution(T min, T max) {
  return uniform_dist<T>(min, max);
}

template <typename T = double>
T exponential(T lambda) {
  // We are not using std::exponential_distribution because of a bug in the
  // implementations by clang and gcc.
  return -std::log(canonical_nonzero()) / lambda;
}

template <typename T = double>
T expo(T A, T x1, T x2) {
  const T a1 = A * x1, a2 = A * x2;
  const T a_min = std::log(std::numeric_limits<T>::min());
#ifndef NDEBUG
  assert(A > T(0.) && x1 >= x2 && a1 > a_min);
#endif
  const T r1 = std::exp(a1);
  const T r2 = a2 > a_min ? std::exp(a2) : T(0.);  // prevent underflow
  T x;
  do {
    /* sample repeatedly until x is in the requested range
     * (it can get outside due to numerical errors, see issue #2959) */
    x = std::log(uniform(r1, r2)) / A;
  } while (!(x <= x1 && x > x2));
  return x;
}

template <typename T>
int sgn(T val) {
  return (T(0) < val) - (val < T(0));
}

template <typename T = double>
T power(T n, T xMin, T xMax) {
  const T n1 = n + 1;
  if (std::abs(n1) < 1E-3) {
    return xMin * std::pow(xMax / xMin, canonical());
  } else if (xMin > 0. && xMax > 0.) {
    return std::pow(uniform(std::pow(xMin, n1), std::pow(xMax, n1)), 1. / n1);
  } else {
    return sgn(xMin) * std::pow(uniform(std::pow(std::abs(xMin), n1),
                                        std::pow(std::abs(xMax), n1)),
                                1. / n1);
  }
}

template <typename T>
int poisson(const T &lam) {
  return std::poisson_distribution<int>(lam)(engine);
}

template <typename T>
int binomial(const int N, const T &p) {
  return std::binomial_distribution<int>(N, p)(engine);
}

template <typename T>
double normal(const T &mean, const T &sigma) {
  return std::normal_distribution<double>(mean, sigma)(engine);
}

template <typename T>
class discrete_dist {
 public:
  discrete_dist() : distribution({1.0}) {}

  explicit discrete_dist(const std::vector<T> &plist)
      : distribution(plist.begin(), plist.end()) {}

  explicit discrete_dist(std::initializer_list<T> l) : distribution(l) {}

  void reset_weights(const std::vector<T> &plist) {
    distribution = std::discrete_distribution<>(plist.begin(), plist.end());
  }
  int operator()() { return distribution(engine); }

 private:
  std::discrete_distribution<> distribution;
};

template <typename T = double>
T cauchy(T pole, T width, T min, T max) {
  /* Use double-precision variables, in order to work around a glibc bug in
   * tanf:
   * https://sourceware.org/bugzilla/show_bug.cgi?id=18221 */
  const double u_min = std::atan((min - pole) / width);
  const double u_max = std::atan((max - pole) / width);
  const double u = uniform(u_min, u_max);
  return pole + width * std::tan(u);
}

template <typename T = double>
T beta(T a, T b) {
  // Otherwise the integral over probability density diverges
  assert(a > T(0.0) && b > T(0.0));
  const T x1 = std::gamma_distribution<T>(a)(engine);
  const T x2 = std::gamma_distribution<T>(b)(engine);
  return x1 / (x1 + x2);
}

template <typename T = double>
T beta_a0(T xmin, T b) {
  assert(xmin > T(0.0) && xmin < T(1.0));
  T y;
  do {
    y = uniform(0.0, -std::log(xmin));
  } while (std::pow((1.0 - std::exp(-y)), b) < canonical());
  return std::exp(-y);
}

class BesselSampler {
 public:
  BesselSampler(const double poisson_mean1, const double poisson_mean2,
                const int fixed_difference);

  std::pair<int, int> sample();

 private:
  static double r_(int n, double a);

  random::discrete_dist<double> dist_;

  double m_;

  const double a_;

  const int N_;

  const bool N_is_positive_;

  static constexpr double m_switch_method_ = 6.0;

  static constexpr double negligible_probability_ = 1.e-12;

  double mu_;

  double sigma_;
};

}  // namespace random
}  // namespace smash

#endif  // SRC_INCLUDE_RANDOM_H_