random.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2012-2020,2026
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
 
#ifndef SRC_INCLUDE_SMASH_RANDOM_H_
#define SRC_INCLUDE_SMASH_RANDOM_H_
 
#include <cassert>
#include <limits>
#include <random>
#include <sstream>
#include <utility>
#include <vector>
 
#include "smash/macros.h"
 
namespace smash {
 
/** Namespace random provides functions for random Number Generation.
 */
 
namespace random {
 
/// The random number engine used is the Mersenne Twister.
using Engine = std::mt19937_64;
 
/// The engine that is used commonly by all distributions.
extern /*thread_local (see commit 897d0b8)*/ Engine engine;
 
/** Provides uniform random numbers on a fixed interval.
 *
 * objects of uniform_dist can be used to provide a large number of
 * random numbers in the same interval. Example:
 *
 * \code
 *   using namespace random;
 *   double sum = 0.0;
 *   auto uniform_0_to_3 = uniform_dist(0., 3.);
 *   for (MANY_TIMES) {
 *     sum += uniform_0_to_3();
 *   }
 * \endcode
 *
 * The random number engine is completely hidden inside the object.
 */
template <typename T>
class uniform_dist {
 public:
  /**
   * Creates the object and fixes the interval.
   *
   * \param min Lower bound of interval.
   * \param max Upper bound of interval.
   * */
  uniform_dist(T min, T max) : distribution(min, max) {}
  /** \returns A random number in the interval. */
  T operator()() { return distribution(engine); }
 
 private:
  /** The distribution object that is being used. */
  std::uniform_real_distribution<T> distribution;
};
 
/** Generates a seed with a truly random 63-bit value, if possible */
int64_t generate_63bit_seed();
 
/** Sets the seed of the random number engine. */
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));
}
 
/// Advance the engine's state and return the generated value.
inline Engine::result_type advance() { return engine(); }
 
/**
 * \returns A uniformly distributed random real number \f$\chi \in [{\rm
 * min}, {\rm max})\f$
 *
 * \param min Minimal sampled value.
 * \param max Maximal sampled value.
 */
template <typename T>
T uniform(T min, T max) {
  /* Strictly speaking, a distribution could be created with min == max, but
   * then it would be UB to use the call operator(), which we do here. */
  assert(min < max);
  return std::uniform_real_distribution<T>(min, max)(engine);
}
 
/**
 * \return A uniformly distributed random integer number \f$\chi \in [{\rm
 * min}, {\rm max})\f$
 *
 * \param min Minimal sampled value.
 * \param max Maximal sampled value.
 */
template <typename T>
T uniform_int(T min, T max) {
  /* Strictly speaking, a distribution could be created with min == max, but
   * then it would be UB to use the call operator(), which we do here. */
  assert(min < max);
  return std::uniform_int_distribution<T>(min, max)(engine);
}
 
/**
 * \return a uniformly distributed random number \f$\chi \in [0,1)\f$.
 *
 * Note that the popular implementations in GCC and clang may return 1:
 *
 * https://gcc.gnu.org/bugzilla/show_bug.cgi?id=64351
 * https://llvm.org/bugs/show_bug.cgi?id=18767
 */
template <typename T = double>
T canonical() {
  return std::generate_canonical<T, std::numeric_limits<double>::digits>(
      engine);
}
 
/**
 * \return A uniformly distributed random number \f$\chi \in (0,1]\f$.
 */
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));
}
 
/**
 * \return A uniform_dist object.
 * \param min Lower bound of interval.
 * \param max Upper bound of interval.
 * */
template <typename T>
uniform_dist<T> make_uniform_distribution(T min, T max) {
  return uniform_dist<T>(min, max);
}
 
/**
 * Draws an exponentially distributed random number.
 *
 * Probability for a given return value \f$\chi\f$ is \f$p(\chi) =
 * \Theta(\chi) \cdot \exp(-t)\f$.
 *
 * \param lambda Rate parameter.
 * \return Sampled random number.
 */
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;
}
 
/**
 * Draws a random number x from an exponential distribution exp(A*x), where A is
 * assumed to be positive, and x is typically negative.
 * The result x is restricted to lie between x1 and x2 (with x2 < x <= x1).
 *
 * \param A Positive shape parameter.
 * \param x1 Maximal sampled value.
 * \param x2 Minimal sampled value.
 * \return Sampled random number.
 *
 * \throw std::logic_error if the computed sampling interval is degenerate,
 *        reversed, or otherwise invalid.
 */
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());
  assert(A > T(0.) && x1 > x2 && a1 > a_min);
  const T high = std::exp(a1);
  const T low = a2 > a_min ? std::exp(a2) : T(0.);  // prevent underflow
  if (unlikely(!(low < high))) {                    // catches NANs, too
    std::ostringstream error_message{};
    error_message << "Function " << __func__
                  << ": internal invariant 'low < high' violated (low = " << low
                  << ", high = " << high << ")";
    throw std::logic_error(error_message.str());
  }
  T x{};
  do {
    /* sample repeatedly until x is in the requested range
     * (it can get outside due to numerical errors). */
    x = std::log(uniform(low, high)) / A;
  } while (!(x <= x1 && x > x2));
  return x;
}
 
/**
 * Signum function.
 *
 * \param val The input value.
 * \return The sign of the input value.
 */
template <typename T>
int sgn(T val) {
  return (T(0) < val) - (val < T(0));
}
 
/**
 
* Sample from a power-law probability density proportional to |x|^n.
* The sample is drawn on the interval [xMin, xMax]. The interval must lie
* entirely on one side of zero; intervals crossing zero are not supported.
* Negative intervals are handled by sampling the absolute values and restoring
* the negative sign.
* For n ≈ -1, the distribution is sampled using the logarithmic limit.
* \tparam T Floating-point type.
* \param n Power-law exponent.
* \param xMin Lower interval bound.
* \param xMax Upper interval bound.
* \return Random value distributed as p(x) ∝ |x|^n on the given interval.
* \throws std::invalid_argument if the interval crosses zero.
* \throws std::invalid_argument if the interval touches zero and n <= -1,
*  where the distribution is not normalizable.
*/
template <typename T = double>
T power(T n, T xMin, T xMax) {
  const T n1 = n + T(1);
 
  if ((xMin < 0 && xMax > 0) || (xMax < 0 && xMin > 0)) {
    throw std::invalid_argument(
        "power: interval crossing zero is not supported");
  } else if ((xMin == 0 || xMax == 0) && n <= -1) {
    throw std::invalid_argument("power: distribution not normalizable at x=0");
  }
 
  if (xMin > xMax) {
    std::swap(xMin, xMax);
  }
  const T sign = xMax < T(0) ? T(-1) : T(1);
 
  const T lo = std::abs(xMin);
  const T hi = std::abs(xMax);
 
  if (std::abs(n1) < T(1e-3)) {
    return sign * lo * std::pow(hi / lo, canonical());
  }
 
  T a = std::pow(lo, n1);
  T b = std::pow(hi, n1);
 
  if (a > b) {
    std::swap(a, b);
  }
 
  return sign * std::pow(uniform(a, b), T(1) / n1);
}
/**
 * Returns a Poisson distributed random number.
 *
 * Probability for a given return value \f$\chi\f$ is \f$p(\chi) =
 * \chi^i/i! \cdot \exp(-\chi)\f$
 *
 * \param lam Mean value of the distribution.
 * \return Sampled random number.
 */
template <typename T>
int poisson(const T &lam) {
  return std::poisson_distribution<int>(lam)(engine);
}
 
/**
 * Returns a binomially distributed random number.
 *
 * \param N Number of trials.
 * \param p Probability of a trial generating true.
 * \return Sampled random number.
 */
template <typename T>
int binomial(const int N, const T &p) {
  return std::binomial_distribution<int>(N, p)(engine);
}
 
/**
 * Returns a random number drawn from a normal distribution.
 *
 * \param mean Mean value of the distribution.
 * \param sigma Standard deviation of the distribution.
 * \return Sampled random number.
 */
template <typename T>
double normal(const T &mean, const T &sigma) {
  return std::normal_distribution<double>(mean, sigma)(engine);
}
 
/**
 * Discrete distribution with weight given by probability vector.
 */
template <typename T>
class discrete_dist {
 public:
  /** Default discrete distribution.
   *
   * Always draws 0.
   */
  discrete_dist() : distribution({1.0}) {}
 
  /** Construct from probability vector.
   * \param plist Vector with probabilities such that P(i) = vec[i]
   */
  explicit discrete_dist(const std::vector<T> &plist)
      : distribution(plist.begin(), plist.end()) {}
 
  /** Construct from probability list.
   * \param l Initializer list with probabilities such that P(i) = l[i]
   */
  explicit discrete_dist(std::initializer_list<T> l) : distribution(l) {}
 
  /** Reset the discrete distribution from a new probability list.
   * \param plist Vector with probabilities such that P(i) = vec[i]
   */
  void reset_weights(const std::vector<T> &plist) {
    distribution = std::discrete_distribution<>(plist.begin(), plist.end());
  }
  /** Draw a random number from the discrete distribution.
   * \return Sampled value
   */
  int operator()() { return distribution(engine); }
 
 private:
  /** The distribution object that is being used. */
  std::discrete_distribution<> distribution;
};
 
/**
 * Draws a random number from a Cauchy distribution (sometimes also called
 * Lorentz or non-relativistic Breit-Wigner distribution) with the given
 * parameters (constant width!) inside the range [min,max]. This function is
 * similar to std::cauchy_distribution, but can return values inside a limited
 * interval.
 * \param pole Pole parameter of the Cauchy function, i.e. location of the peak.
 * \param width Width parameter of the Cauchy function, determining the
 * sharpness of the peak.
 * \param min Minimum value to be returned.
 * \param max Maximum value to be returned.
 * \return Sampled random number.
 */
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);
}
 
/**
 * Draws a random number from a beta-distribution, where probability density of
 * \f$x\f$ is \f$p(x) = frac{\Gamma(a)\Gamma(b)}{Gamma(a+b)}
 *  x^{a-1} (1-x)^{b-1}\f$. This distribution is necessary for string
 *  formation. The implementation uses a property connecting beta distribution
 * to gamma-distribution. Interchanging a and b will not change results.
 *
 * \param a Shape parameter.
 * \param b Scale parameter.
 * \return Sampled random number.
 */
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);
}
 
/**
 * Draws a random number from a beta-distribution with a = 0. In this case
 * the probability density is \f$p(x) = 1/x (1-x)^b\f$. The integral from
 * 0 to 1 over this distribution diverges, so the sampling is performed
 * in the interval (xmin, 1). This distribution is necessary for string
 * formation. The implementation uses the following property:
 * \f$p(x)dx = dx/x (1-x)^b = (1-x)^b d ln(x) = (1 - e^{-y})^b dy\f$, where
 * \f$ y = - ln(x) \f$.
 *
 * \param xmin Minimal sampled value.
 * \param b Second shape parameter.
 * \return Sampled random number.
 */
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);
}
 
/**
 * The intention of this class is to efficiently sample \f$ (N_1, N_2) \f$
 * from the Bessel distribution \f$ p(N_1,N_2) \sim \mathrm{Poi}(\nu_1)
 * \mathrm{Poi}(\nu_2) \delta(N_1 - N_2 = N)\f$, where \f$\mathrm{Poi}(\nu)\f$
 * denotes Poisson distribution with mean \f$\nu\f$. In other words, this
 * class samples two Poisson numbers with a given mean and a fixed difference.
 * The intended use is to sample the number of baryons and antibaryons, given
 * their means and net baryon number.
 *
 * The distribution of \f$ min(N_1,N_2) \f$ is a so-called Bessel distribution.
 * Denoting \f$ a = \sqrt{\nu_1 \nu_2}\f$, \f$ p(N_{smaller} = k) =
 * \frac{(a/2)^{2k+N}}{I_N(a) k! (N+k)!} \f$. We sample this distribution using
 * the method suggested by Yuan and Kalbfleisch \cite Yuan2000 : if
 * \f$ m = \frac{1}{2} (\sqrt{a^2 + N^2} - N) > 6\f$, then the distribution is
 * approximated well by a Gaussian, else probabilities are computed explicitely
 * and a table sampling is applied.
 */
class BesselSampler {
 public:
  /**
   * Construct a \ref BesselSampler.
   *
   * \param[in] poisson_mean1 Mean of the first number's Poisson distribution.
   * \param[in] poisson_mean2 Mean of the second number's Poisson distribution.
   * \param[in] fixed_difference Difference between the sampled numbers.
   * \return Constructed sampler.
   */
  BesselSampler(const double poisson_mean1, const double poisson_mean2,
                const int fixed_difference);
 
  /**
   * Sample two numbers from given Poissonians with a fixed difference.
   *
   * \return Pair of first and second sampled number.
   */
  std::pair<int, int> sample();
 
 private:
  /**
   * Compute the ratio of two Bessel functions
   * r(n,a) = bessel_I(n+1,a)/bessel_I(n,a) using the continued fraction
   * representation (see \cite Yuan2000).
   *
   * \param[in] n First Bessel parameter.
   * \param[in] a Second Bessel parameter.
   * \return Ratio bessel_I(n+1,a)/bessel_I(n,a).
   */
  static double r_(int n, double a);
 
  /// Vector to store tabulated values of probabilities for small m case (m <6).
  random::discrete_dist<double> dist_;
 
  /// Mode of the Bessel function, see \cite Yuan2000 for details.
  double m_;
 
  /// Second parameter of Bessel distribution, see \cite Yuan2000 for details.
  const double a_;
 
  /// First parameter of Bessel distribution (= \f$ \nu \f$ in \cite Yuan2000).
  const int N_;
 
  /// Boolean variable to verify that N > 0.
  const bool N_is_positive_;
 
  /**
   * Switching mode to normal approximation.
   * \note Normal approximation of Bessel functions is possible for modes >= 6.
   * See \cite Yuan2000 for details.
   */
  static constexpr double m_switch_method_ = 6.0;
 
  /// Probabilities smaller than negligibly_probability are neglected.
  static constexpr double negligible_probability_ = 1.e-12;
 
  /// Mean of the Bessel distribution.
  double mu_;
 
  /// Standard deviation of the Bessel distribution.
  double sigma_;
};
 
}  // namespace random
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_RANDOM_H_