Random number generation¶

The 2011 ISO C++ standard defined a new functional paradigm for sampling from random number distributions. In this paradigm, random number engines generate a uniformly distributed stream of bits. Then, distributions use that entropy to sample a random number from a distribution.

Engines¶

Celeritas defaults to using an in-house implementation of the XORWOW [Marsaglia, 2003] bit shifting generator. Each thread’s state is seeded at runtime by filling the state with bits generated from a 32-bit Mersenne twister. When a new event begins through the Geant4 interface, each thread’s state is initialized using same seed and skipped ahead a different number of subsequences so the sequences on different threads will not have statistically correlated values.

void celeritas::initialize_xorwow(Span<XorwowState> state, XorwowSeed const &seed, StreamId stream)¶

class XorwowRngEngine¶

Generate random data using the XORWOW algorithm.

The XorwowRngEngine uses a C++11-like interface to generate random data. The sampling of uniform floating point data is done with specializations to the GenerateCanonical class.

The resize function for XorwowRngStateData will fully randomize the state at initialization. Alternatively, the state can be initialized with a seed, subsequence, and offset.

See Marsaglia (2003) for the theory underlying the algorithm and the the “example” xorwow that combines an xorshift output with a Weyl sequence (https://www.jstatsoft.org/index.php/jss/article/view/v008i14/916).

Initialization moves the state ahead to the given subsequence (a subsequence has size \(2^{67}\)) and skips offset random numbers. It is recommended to initialize the state using a very different generator from the one being initialized to avoid correlations. Here, the 64-bit SplitMix64 generator is used for initialization [Matsumoto et al., 2007] .

For a description of the jump ahead method using the polynomial representation of the recurrence, see Haramoto et al. [2008] .

The jump polynomials were precomputed using a python script. For a more detailed description of how to calculate the jump polynomials using Knuth’s square-and-multiply algorithm in \( O(k^2 log d) \) time (where k is the number of bits in the state and d is the number of steps to skip ahead), see Collins [2008] .

Distributions¶

Distributions are function-like objects whose constructors take the parameters of the distribution: for example, a uniform distribution over the range \([a, b)\) takes the a and b parameters as constructor arguments. The templated call operator accepts a random engine as its sole argument.

Celeritas extends this paradigm to physics distributions. At a low level, it has random number distributions that result in single real values (such as uniform, exponential, gamma) and correlated three-vectors (such as sampling an isotropic direction).

class BernoulliDistribution¶

Select one of two options with a given probability.

The constructor argument is the chance of returning true, with an optional second argument for normalizing with a fraction of false values.

BernoulliDistribution snake_eyes(1, 35);
if (snake_eyes(rng))
{
    // ...
}
BernoulliDistribution also_snake_eyes(1.0 / 36.0);

template<class T> class DeltaDistribution¶: Distribution where the “sampled” value is just the given value.

template<class RealType = ::celeritas::real_type> class ExponentialDistribution¶

Sample from an exponential distribution.

Sample from a probability distribution function with the normalized PDF:

\[ f(x; \lambda) = \lambda e^{-\lambda x} \quad \mathrm{for} \ x >= 0 \]

which integrated into a CDF and inverted gives a sample:

\[ x = -\frac{\log \xi}{ \lambda} \]

This is simply an implementation of std::exponential_distribution that uses the Celeritas canonical generator and is independent of library implementation.

Note (for performance-critical sections of code) that if this class is constructed locally with the default value of lambda = 1.0, the inversion and multiplication will be optimized out.

template<class RealType = ::celeritas::real_type> class GammaDistribution¶

Sample from a gamma distribution.

The gamma distribution can be parameterized with a shape parameter \( \alpha \) and a scale parameter \( \beta \) and has the PDF:

\[ f(x; \alpha, \beta) = \frac{x^{\alpha - 1} e^{-x / \beta}}{\beta^\alpha \Gamma(\alpha)} \quad \mathrm{for}\ x > 0, \quad \alpha, \beta > 0 \]

The algorithm described in Marsaglia and Tsang [2000] is used here to generate gamma-distributed random variables. The steps are:

Set \( d = \alpha - 1/3 \) and \( c = 1 / \sqrt{9d} \)
Generate random variates \( Z \sim N(0,1) \) and \( U \sim U(0,1) \)
Set \( v = (1 + cZ)^3 \)
If \( \log U < Z^2 / 2 + d(1 - v + \log v) \) return \( dv \). Otherwise, go to step 2. A squeeze function can be used to avoid the two logarithms in most cases by accepting early if \( U < 1 - 0.0331 Z^4 \).

Though this method is valid for \( \alpha \ge 1 \), it can easily be extended for \( \alpha < 1 \): if \( X \sim \Gamma(\alpha + 1) \) and \( U \sim U(0,1) \), then \( X U^{1/\alpha} \sim \Gamma(\alpha) \).

template<class RealType = ::celeritas::real_type> class InverseSquareDistribution¶

Sample \( 1/x^2 \) over a given domain.

This distribution is defined on a positive range \( [a, b) \) and has the normalized PDF:

\[ f(x; a, b) = \frac{a b}{x^2 (b - a)} \quad \mathrm{for} \ a \le x < b \]

which integrated into a CDF and inverted gives a sample:

\[ x = \frac{a b}{a + \xi (b - a)} \]

template<class RealType = ::celeritas::real_type> class IsotropicDistribution¶: Sample points uniformly on the surface of a unit sphere.

template<class RealType = ::celeritas::real_type> class NormalDistribution¶

Sample from a normal distribution.

The normal (Gaussian) distribution is

\[ p(x) = \exp\!\left( -\frac{(x - \mu)^{2}}{2\sigma^{2}} \right) \]

with mean \(\mu\) and standard deviation \(\sigma\).

This implementation uses the Box-Muller transform to generate pairs of independent, normally distributed random numbers, and returns them one at a time. Two random numbers uniformly distributed on [0, 1] are mapped to two independent, standard, normally distributed samples using the relations:

\[ x_1 = \sqrt{-2 \ln \xi_1} \cos(2 \pi \xi_2) x_2 = \sqrt{-2 \ln \xi_1} \sin(2 \pi \xi_2) \]

template<class RealType = ::celeritas::real_type> class PoissonDistribution¶

Sample from a Poisson distribution.

The Poisson distribution describes the probability of \( k \) events occurring in a fixed interval given a mean rate of occurrence \( \lambda \) and has the PMF:

\[ f(k; \lambda) = \frac{\lambda^k e^{-\lambda}}{k!} \:. \]

For small \( \lambda \), a direct method described in [Knuth, 1968] can be used to generate samples from the Poisson distribution. Uniformly distributed random numbers are generated until the relation

\[ \prod_{k = 1}^n U_k \le e^{-\lambda} \]

is satisfied; then, the random variable \( X = n - 1 \) will have a Poisson distribution. On average this approach requires the generation of \( \lambda + 1 \) uniform random samples, so a different method should be used for large \( \lambda \).

Geant4 uses Knuth’s algorithm for \( \lambda \le 16 \) and a Gaussian approximation for \( \lambda > 16 \) (see G4Poisson), which is faster but less accurate than other methods. The same approach is used here.

Todo:: Break this into two distributions: one actual poisson distribution, one “integer normal” distribution, and a variant type that selects between them. In most cases we care about, lambda is small.

template<class RealType = ::celeritas::real_type> class PowerDistribution¶

Sample a power distribution for powers not equal to -1.

This distribution for a power \( p != -1 \) is defined on a positive range \( [a, b) \) and has the normalized PDF:

\[ f(x; p, a, b) = x^{p}\frac{p + 1}{b^{p + 1} - a^{p + 1}} \quad \mathrm{for } a < x < b \]

which integrated into a CDF and inverted gives a sample:

\[ x = \left[ (b^{p + 1} - a^{p + 1})\xi + a^{p + 1} \right]^{1/(p + 1)} \]

The quantity in brackets is a uniform sample on \( [a^{p + 1}, b^{p + 1}) \)

See C13A (and C15A) from Everett and Cashwell [1972] , with \( x \gets u \), \( p \gets m - 1 \).

Note

For \( p = -1 \) see ReciprocalDistribution , and in the degenerate case of \( p = 0 \) this is mathematically equivalent to UniformRealDistribution .

template<class RealType = ::celeritas::real_type> class RadialDistribution¶: Sample from a uniform radial distribution.

template<class RealType = ::celeritas::real_type> class ReciprocalDistribution¶

Reciprocal or log-uniform distribution.

This distribution is defined on a positive range \( [a, b) \) and has the normalized PDF:

\[ f(x; a, b) = \frac{1}{x (\ln b - \ln a)} \quad \mathrm{for} \ a \le x < b \]

which integrated into a CDF and inverted gives a sample:

\[ x = a \left( \frac{b}{a} \right)^{\xi} = a \exp\!\left(\xi \log \frac{b}{a} \right) \]

template<class RealType = ::celeritas::real_type> class UniformBoxDistribution¶: Sample a point uniformly in a box.

template<class RealType = ::celeritas::real_type> class UniformRealDistribution¶

Sample from a uniform distribution.

This distribution is defined between two arbitrary real numbers a and b , and has a flat PDF between the two values. It is allowable for the two numbers to have reversed order. The normalized PDF is:

\[ f(x; a, b) = \frac{1}{b - a} \quad \mathrm{for} \ a \le x < b \]

which integrated into a CDF and inverted gives a sample:

\[ x = (b - a) \xi + a \]

Additionally we define a few helper classes for common physics sampling routines.

template<class RealType = ::celeritas::real_type> class RejectionSampler¶

Return whether a rejection loop needs to continue trying.

A common implementation of sampling from a “difficult” (non-analytically invertible) probability distribution function is to bound the difficult distribution f(x) with another easily sampled function g(x) . Given a maximum value M over the x interval being sampled, it is equivalent to sampling f(x) by instead sampling from g(x) and rejecting with probability

\[ \frac{f(x)}{M g(x)} \]

These invocations generate statistically equivalent results:

BernoulliDistribution(1 - p / pmax)(rng);
!BernoulliDistribution(p / pmax)(rng);
p < pmax * UniformDistribution{}(rng)
RejectionSampler(p, pmax)(rng);

This is meant for rejection sampling, e.g., on cross section:

do {
  xs = sample_xs(rng);
} while (RejectionSampler{xs, xs_max}(rng));

template<class F, class T> class Selector¶

On-the-fly selection of a weighted discrete distribution.

This algorithm encapsulates the loop for sampling from distributions of a function f(index) -> real (usually with the function prototype real_type(*)(size_type) ). The index can be an integer, an enumeration, or an OpaqueId . The Selector is constructed with the size of the distribution (but using the indexing type).

Edge cases are thoroughly tested (it will never iterate off the end if normalized, even for incorrect values of the “total” probability/xs), and it uses one fewer register than the typical accumulation algorithm. When building with debug checking and normalization, the constructor asserts that the provided “total” value is consistent.

The given function must return a consistent value for the same given argument.

auto select_el = make_selector(
    [](ElementId i) { return xs[i.get()]; },
    ElementId{num_elements()},
    tot_xs);
ElementId el = select_el(rng);

auto select_val = make_selector([](size_type i) { return pdf[i]; },
                                pdf.size());
size_type idx = select_val(rng);

Create a normalized selector from a function and total accumulated value.