<jupyter_start><jupyter_text>PyData Amsterdam 2019 :: Numba Tutorial Notebook 6 :: Extra Topics Random NumbersGPUs can be extremely useful for Monte Carlo applications where you need to use large amounts of random numbers.  CUDA ships with an excellent set of random number generation algorithms in the cuRAND library.  Unfortunately, cuRAND is defined in a set of C headers which Numba can't easily compile or link to.  (Numba's CUDA JIT does not ever create C code for CUDA kernels.)  It is on the Numba roadmap to find a solution to this problem, but it may take some time.In the meantime, Numba version 0.33 and later includes the `xoroshiro128+` generator, which is pretty high quality, though with a smaller period ($2^{128} - 1$) than the XORWOW generator in cuRAND.  To use it, you will want to initialize the RNG state on the host for each thread in your kernel:<jupyter_code>import numpy as np
from numba import cuda
from numba.cuda.random import create_xoroshiro128p_states, xoroshiro128p_uniform_float32

threads_per_block = 64
blocks = 24
rng_states = create_xoroshiro128p_states(threads_per_block * blocks, seed=1)<jupyter_output><empty_output><jupyter_text>This state creation function initializes each state to be in the same sequence designated by the seed, but separated by $2^{64}$ steps from each other.  This ensures that different threads will not accidentally end up with overlapping sequences (since you will not have enough patience to draw anywhere near $2^{64}$ random numbers in a single thread).We can use these random number states in our kernel by passing it in as an argument:<jupyter_code>@cuda.jit
def compute_pi(rng_states, iterations, out):
    """Find the maximum value in values and store in result[0]"""
    thread_id = cuda.grid(1)

    # Compute pi by drawing random (x, y) points and finding what
    # fraction lie inside a unit circle
    inside = 0
    for i in range(iterations):
        x = xoroshiro128p_uniform_float32(rng_states, thread_id)
        y = xoroshiro128p_uniform_float32(rng_states, thread_id)
        if x**2 + y**2 <= 1.0:
            inside += 1

    out[thread_id] = 4.0 * inside / iterations
out = np.zeros(threads_per_block * blocks, dtype=np.float32)
compute_pi[blocks, threads_per_block](rng_states, 10000, out)
print('pi:', out.mean())<jupyter_output><empty_output><jupyter_text>Shared MemoryWe briefly mention in notebook 4 that the CUDA programming model organizes threads into a two-layer structure.  A grid is composed of many blocks, which are composed of many threads.  Threads within the same block can communicate much more easily than threads in different blocks.  The main mechanism for this communication is *shared memory*.  Shared memory is discussed extensively in the CUDA C Programming Guide, as well as many other books on CUDA programming.  We will only describe it very briefly here, and focus mainly on the Python syntax for using it.Shared memory is a section of memory that is visible at the block level.  Different blocks cannot see each other's shared memory, and all the threads within a block see the same shared memory.  It does not persist after a CUDA kernel finishes executing.  Shared memory is scarce hardware resource, so should be used sparingly or side effects such as lower performance or even kernel launch failure (if you exceed the hardware limit of 48 kB per block) will occur.Shared memory is good for several things: * caching of lookup tables that will be randomly accessed * buffering output from threads so it can be coalesced before writing it back to device memory. * staging data for scatter/gather operations within a block As an example of the power of shared memory, let's write a transpose kernel that takes a 2D array in row-major order and puts it in column-major order.  (This is based on Mark Harris' blog post at: https://devblogs.nvidia.com/parallelforall/efficient-matrix-transpose-cuda-cc/)First, let's do the naive approach where we let each thread read and write individual elements independently:<jupyter_code>TILE_DIM = 32
BLOCK_ROWS = 8

@cuda.jit
def transpose(a_in, a_out):
    x = cuda.blockIdx.x * TILE_DIM + cuda.threadIdx.x
    y = cuda.blockIdx.y * TILE_DIM + cuda.threadIdx.y

    for j in range(0, TILE_DIM, BLOCK_ROWS):
        a_out[x, y + j] = a_in[y + j, x]
size = 1024
a_in = cuda.to_device(np.arange(size*size, dtype=np.int32).reshape((size, size)))
a_out = cuda.device_array_like(a_in)

print(a_in.copy_to_host())
grid_shape = (int(size/TILE_DIM), int(size/TILE_DIM))
%timeit transpose[grid_shape,(TILE_DIM, BLOCK_ROWS)](a_in, a_out); cuda.synchronize()
print(a_out.copy_to_host())<jupyter_output><empty_output><jupyter_text>Now let's use shared memory to copy a 32x32 tile at a time.  We'll use a global value for the tile size so it will be known act compile time:<jupyter_code>import numba.types

TILE_DIM_PADDED = TILE_DIM + 1  # Read <NAME>' blog post to find out why this improves performance!

@cuda.jit
def tile_transpose(a_in, a_out):
    # THIS CODE ASSUMES IT IS RUNNING WITH A BLOCK DIMENSION OF (TILE_SIZE x TILE_SIZE)
    # AND INPUT IS A MULTIPLE OF TILE_SIZE DIMENSIONS
    tile = cuda.shared.array((TILE_DIM, TILE_DIM_PADDED), numba.types.int32)

    x = cuda.blockIdx.x * TILE_DIM + cuda.threadIdx.x
    y = cuda.blockIdx.y * TILE_DIM + cuda.threadIdx.y
    
    for j in range(0, TILE_DIM, BLOCK_ROWS):
        tile[cuda.threadIdx.y + j, cuda.threadIdx.x] = a_in[y + j, x] # transpose tile into shared memory

    cuda.syncthreads()  # wait for all threads in the block to finish updating shared memory

    #Compute transposed offsets
    x = cuda.blockIdx.y * TILE_DIM + cuda.threadIdx.x
    y = cuda.blockIdx.x * TILE_DIM + cuda.threadIdx.y

    for j in range(0, TILE_DIM, BLOCK_ROWS):
        a_out[y + j, x] = tile[cuda.threadIdx.x, cuda.threadIdx.y + j];

a_out = cuda.device_array_like(a_in) # replace with new array

%timeit tile_transpose[grid_shape,(TILE_DIM, BLOCK_ROWS)](a_in, a_out); cuda.synchronize()
print(a_out.copy_to_host())<jupyter_output><empty_output><jupyter_text>That's a 40% speed up!  Generalized UfuncsUfuncs broadcast a scalar function over array inputs but what if you want to broadcast a lower dimensional array function over a higher dimensional array?  This is called a *generalized ufunc* ("gufunc"), and it opens up a whole new frontier for applying ufuncs.Generalized ufuncs are a little more tricky becuase they need a *signature* (not to be confused with the Numba type signature) that shows the index ordering when dealing with multiple inputs.  Fully explaining "gufunc" signatures is beyond the scope of this tutorial, but you can learn more from:* The NumPy docs on gufuncs: https://docs.scipy.org/doc/numpy/reference/c-api.generalized-ufuncs.html* The Numba docs on gufuncs: http://numba.pydata.org/numba-doc/latest/user/vectorize.htmlthe-guvectorize-decorator* The Numba docs on CUDA gufuncs: http://numba.pydata.org/numba-doc/latest/cuda/ufunc.htmlgeneralized-cuda-ufuncsLet's write our own normalization function.  This will take an array input and compute the L2 norm along the last dimension.  Generalized ufuncs take their output array as the last argument, rather than returning a value.  If the output is a scalar (as it will be for the norm function), then we will still receive a 1D array as the output argument, but we will write the scalar output to element 0 of the array.<jupyter_code>from numba import guvectorize
import math

@guvectorize(['(float32[:], float32[:])'], # have to include the output array in the type signature
             '(i)->()',                 # map a 1D array to a scalar output
             target='cuda')
def l2_norm(vec, out):
    acc = 0.0
    for value in vec:
        acc += value**2
    out[0] = math.sqrt(acc)<jupyter_output><empty_output><jupyter_text>To test this, let's construct some points on the unit circle:<jupyter_code>angles = np.random.uniform(-np.pi, np.pi, 10)
coords = np.stack([np.cos(angles), np.sin(angles)], axis=1)
print(coords)<jupyter_output>[[ 0.87448536 -0.48505191]
 [-0.65614478  0.75463503]
 [-0.12792447 -0.99178391]
 [ 0.85177927  0.52390082]
 [ 0.58268136 -0.81270071]
 [ 0.94916918 -0.31476637]
 [-0.90277458  0.43011401]
 [ 0.13394428  0.99098886]
 [ 0.3989953   0.91695297]
 [ 0.20003634  0.97978848]]<jupyter_text>As expected, the L2 norm is 1.0, up to rounding errors:<jupyter_code>l2_norm(coords)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07990976015198291, 'fraction_numerical': 0.03431488957492282, 'mean_word_length': 4.406290115532734, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 30, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '248076', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 2890, 'n_tokens_neox': 2653, 'n_words': 1031, 'path': '6 - Extra Topics.ipynb', 'repo_name': 'alanlujan91/pydata-amsterdam2019-numba', 'size': '8427'}