Random Number Generation

Random Number Generation#

Pystencils offers some limited pseudo-random number generation facilities for use in generated kernels. These reside in the module pystencils.sympyextensions.random and can be imported through pystencils.random. At the moment, pystencils implements the Philox [1] counter-based random number generator. It is capable of producing independent streams of random numbers in parallel and is therefore feasible for use in parallel CPU and GPU kernels.

Calling RNGs in a kernel#

To use the random number generator in a kernel, first create an instance of ps.random.Philox:

rng = ps.random.Philox("rng", "float32", seed=0xCAFE)

You must pass in a name that is a valid C identifier, and a floating-point data type. The random numbers generated will be of that data type. Furthermore, a 32-bit unsigned integer seed is required.

Use the get_random_vector method to get a symbolic invocation of the generator. Each invocation of the Philox RNG returns a vector of random numbers: four float32, or two float64 values. The method get_random_vector must be called with a counter argument which can be used to distinguish different kernel calls at runtime (e.g. to generate different random numbers at each timestep of a transient numerical solver). It returns two values: an indexable symbol (as an instance of sp.IndexedBase) that represents the random vector, and an assignment which sets that symbol to the result of an RNG invocation. The assignment must be inserted into your kernel at an appropriate place. To access the random numbers, use the []-operator on the symbol. The following example illustrates this by adding up all four generated float32-values and writing them to a field:

t = ps.TypedSymbol("t", "uint32")
rx, rasm = rng.get_random_vector(t)

f = ps.fields("f: float32[2D]")

asms = [
    rasm,
    ps.Assignment(f(), rx[0] + rx[1] + rx[2] + rx[3])
]

asms
\[\displaystyle \left[ rng_{0} \leftarrow_{} \operatorname{rng}_{0}{\left(t \right)}, \ {f}_{(0,0)} \leftarrow_{} {rng_{0}}_{0} + {rng_{0}}_{1} + {rng_{0}}_{2} + {rng_{0}}_{3}\right]\]

We can now generate code for and execute the kernel:

ker = ps.create_kernel(asms)

kfunc = ker.compile()
arr = np.zeros((12, 31), dtype=np.float32)
kfunc(f=arr, t=12)

plt.imshow(arr)

<matplotlib.image.AxesImage at 0x6ffd9274f400>
../_images/14f07d14da777bde56b45217837eca54ddf78b0f3e929e9d174d2bad5171ed2c.png

Ta-da! Random noise.

Controlling RNG Behavior#

The random numbers generated by the Philox RNG are fully determined by six unsigned 32-bit integer numbers. Two of them can be directly user-specified:

  • The seed, which may be specified when creating the Philox instance; the default value is zero.

  • The timestep counter, which must be passed to get_random_vector.

The seed, together with an internal invocation counter (incremented by one at each call to get_random_vector), make up the RNG keys. The timestep, plus the x-, y-, and z-indices of the current kernel iteration, combine to make the four counters. Due to the incrementing hidden key, multiple invocations of the same RNG within a kernel will produce different random numbers. Invoking two different Philox RNGs with the same seed in the same kernel, on the other hand, will result in two identical streams of random numbers.

Contents