Random Numbers
==============

Parallel codes often need random numbers for Monte Carlo methods, randomized initialization, sampling, or synthetic test data.
*alpaka* provides random engines and distributions that can be used directly inside kernels.

Using Random Numbers in a Kernel
--------------------------------

To avoid correlations between random numbers generated by different threads each thread should get its own deterministic engine state generated by a unique seed.
The easiest way to do that is to derive the seed from the loop index.
The unique engine state can then be used to generate as many random samples as needed from the chosen distribution.

*alpaka* supports the following distributions:

- ``UniformReal<float or double>`` for samples in a bounded floating-point interval where all values can occur with equal probability.
- ``NormalReal<float or double>`` for Gaussian-distributed samples with a chosen mean and standard deviation.

Uniform Random Numbers
``````````````````````

  .. literalinclude:: ../../snippets/example/170_random.cpp
    :language: cpp
    :start-after: BEGIN-TUTORIAL-randomKernel
    :end-before: END-TUTORIAL-randomKernel
    :dedent:

This example the engine ``rand::engine::Philox4x32x10`` and used the distribution ``rand::distribution::UniformReal<float>`` with the half-open interval ``[0, 1)`` interval ``rand::interval::co``.
The configuration of intervals is explained in the :ref:`random_number_intervals` section.
The distribution is used to create a random index per element in the output data container.

Launching the Kernel:

  .. literalinclude:: ../../snippets/example/170_random.cpp
    :language: cpp
    :start-after: BEGIN-TUTORIAL-randomLaunch
    :end-before: END-TUTORIAL-randomLaunch
    :dedent:


Uniform distributions are mostly used for probabilities, random offsets, randomized initialization, and rejection sampling.
Normal distributions are used for noise models, perturbations around a mean value, and many Monte Carlo methods.

Normal Distribution
```````````````````

``NormalReal`` generates Gaussian noise with a chosen mean and standard deviation.
Unlike the uniform distribution, it keeps internal state, therefore you should not share the distribution objects between threads.

  .. literalinclude:: ../../snippets/example/170_random.cpp
    :language: cpp
    :start-after: BEGIN-TUTORIAL-randomNormalKernel
    :end-before: END-TUTORIAL-randomNormalKernel
    :dedent:

Launching the kernel is the same as before; only the kernel logic changes.

  .. literalinclude:: ../../snippets/example/170_random.cpp
    :language: cpp
    :start-after: BEGIN-TUTORIAL-randomNormalLaunch
    :end-before: END-TUTORIAL-randomNormalLaunch
    :dedent:

.. _random_number_intervals:

Intervals
---------

``UniformReal`` supports four interval tags:

- ``rand::interval::co`` gives ``[a, b)``
- ``rand::interval::oc`` gives ``(a, b]``
- ``rand::interval::cc`` gives ``[a, b]``
- ``rand::interval::oo`` gives ``(a, b)``

The following kernel shows all four forms side by side.

  .. literalinclude:: ../../snippets/example/170_random.cpp
    :language: cpp
    :start-after: BEGIN-TUTORIAL-randomIntervalsKernel
    :end-before: END-TUTORIAL-randomIntervalsKernel
    :dedent:


Monte Carlo Pi
--------------

A classic example is `Monte Carlo <https://en.wikipedia.org/wiki/Monte_Carlo_method>`__ estimation of pi.
Draw points in the square ``[0, 1) x [0, 1)``, count how many land inside the unit quarter circle, and estimate ``pi`` from that ratio.
The half-open interval matches array-style data access avoids awkward endpoint corner cases.

In the example each worker draws one point, writes ``1`` if the point falls inside the quarter circle, and then a reduction adds up all hits.

  .. literalinclude:: ../../snippets/example/170_random.cpp
    :language: cpp
    :start-after: BEGIN-TUTORIAL-piKernel
    :end-before: END-TUTORIAL-piKernel
    :dedent:

The reduction happens on the same queue right after the kernel.

  .. literalinclude:: ../../snippets/example/170_random.cpp
    :language: cpp
    :start-after: BEGIN-TUTORIAL-piLaunch
    :end-before: END-TUTORIAL-piLaunch
    :dedent:

After copying back the single reduction result, the estimate itself is just the usual Monte Carlo formula.

  .. literalinclude:: ../../snippets/example/170_random.cpp
    :language: cpp
    :start-after: BEGIN-TUTORIAL-piEstimate
    :end-before: END-TUTORIAL-piEstimate
    :dedent:


Complete Source File
--------------------

.. raw:: html

   <details class="full-source">
   <summary>170_random.cpp</summary>

.. filteredliteralinclude:: ../../snippets/example/170_random.cpp
   :language: cpp
   :linenos:

.. raw:: html

   </details>
   <br/>