---
jupytext:
  formats: md:myst
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.4
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
mystnb:
  execution_mode: cache
---

```{code-cell} ipython3
:tags: [remove-cell]

import sympy as sp
import pystencils as ps
import numpy as np
import matplotlib.pyplot as plt

sp.init_printing(use_latex="mathjax")
```

(guide_random_numbers)=
# 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* [^philox] 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 {any}`ps.random.Philox <pystencils.sympyextensions.random.Philox>`:

```{code-cell} ipython3
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 {any}`sp.IndexedBase <sympy.tensor.indexed.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:

```{code-cell} ipython3
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
```

We can now generate code for and execute the kernel:

```{code-cell} ipython3
ker = ps.create_kernel(asms)

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

plt.imshow(arr)
```

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.

[^philox]: J. K. Salmon, M. A. Moraes, R. O. Dror and D. E. Shaw, 
    "Parallel random numbers: As easy as 1, 2, 3,"
    SC '11: Proceedings of 2011 International Conference for High Performance Computing,
    Networking, Storage and Analysis, Seattle, WA, USA, 2011, pp. 1-12,
    doi: [10.1145/2063384.2063405](https://doi.org/10.1145/2063384.2063405).