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 [Mar03] 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)

Initialize XORWOW states with well-distributed random data.

cuRAND and hipRAND implement functions that permute the original (but seemingly arbitrary) seeds given in Marsaglia with 64 bits of seed data. We simply generate pseudorandom, independent starting states for all data in all threads using a 32-bit MT19937 engine.

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).

For a description of the jump ahead method using the polynomial representation of the recurrence, see: Haramoto, H., Matsumoto, M., Nishimura, T., Panneton, F., L’Ecuyer, P. 2008. “Efficient jump ahead for

F2-linear random number generators”. INFORMS Journal on Computing.

https://pubsonline.informs.org/doi/10.1287/ijoc.1070.0251.

The jump polynomials were precomputed using https://github.com/celeritas-project/utils/blob/main/prng/xorwow-jump.py. 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, J. 2008. “Testing, Selection, and Implementation of Random

Number Generators”. ARL-TR-4498.

https://apps.dtic.mil/sti/pdfs/ADA486637.pdf.

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 [MT00] is used here to generate gamma-distributed random variables. The steps are:

1. Set \( d = \alpha - 1/3 \) and \( c = 1 / \sqrt{9d} \)

2. Generate random variates \( Z \sim N(0,1) \) and \( U \sim U(0,1) \)

3. Set \( v = (1 + cZ)^3 \)

4. 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.

This uses the Box-Muller transform to generate pairs of independent, normally distributed random numbers, returning 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 occuring in a fixed interval given a mean rate of occurance \( \lambda \) and has the PMF:

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

For small \( \lambda \), a direct method described in Knuth, Donald E., Seminumerical Algorithms, The Art of Computer Programming, Volume 2 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.

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));

class ElementSelector

Make a weighted random selection of an element.

The element chooser is for selecting an elemental component (atom) of a material based on the microscopic cross section and the abundance fraction of the element in the material.

On construction, the element chooser uses the provided arguments to precalculate all the microscopic cross sections in the given storage space. The given function calc_micro_xs must accept a ElementId and return a real_type, a non-negative microscopic cross section.

The element chooser does not calculate macroscopic cross sections because they’re multiplied by fraction, not number density, and we only care about the fractional abundances and cross section weighting.

ElementSelector select_element(mat, calc_micro, storage);
real_type total_macro_xs
    = select_element.material_micro_xs() * mat.number_density();
ElementComponentId id = select_element(rng);
real_type selected_micro_xs
    = select_element.elemental_micro_xs()[el.get()];
ElementView el = mat.make_element_view(id);
// use el.Z(), etc.

Note that the units of the calculated microscopic cross section will be identical to the units returned by calc_micro_xs. The macroscopic cross section units (micro times mat.number_density() ) will be 1/len if and only if calc_micro units are len^2.

Todo:

Refactor to use Selector.

class IsotopeSelector

Make a weighted random selection of an isotope.

class TabulatedElementSelector

Make a weighted random selection of an element.

This selects an elemental component (atom) of a material based on the precalculated cross section CDF tables of the elements in the material. Unlike ElementSelector which calculates the microscopic cross sections on the fly, this interpolates the values using tabulated CDF grids.

The physics model implementations are built on top of these helper distributions.